tidyverseとNSE
tidyverse系で変数を使いたいときめんどくさいなぁという話。
はじめに要点をざっくりまとめると、tidyverse系でdplyrでdata$列名ではなく列名と直接書けるのはNSE表現になっているので、そういう場所はNSEに合わせた書き方しようね。変数はSEなのでそのままでは評価されないよ、という話。
以下で考える。
library(tidyverse) df = as_tibble(iris) df %>% filter(Species == 'versicolor') %>% head() # # A tibble: 6 x 5 # Sepal.Length Sepal.Width Petal.Length Petal.Width Species # <dbl> <dbl> <dbl> <dbl> <fct> # 1 7 3.2 4.7 1.4 versicolor # 2 6.4 3.2 4.5 1.5 versicolor # 3 6.9 3.1 4.9 1.5 versicolor # 4 5.5 2.3 4 1.3 versicolor # 5 6.5 2.8 4.6 1.5 versicolor # 6 5.7 2.8 4.5 1.3 versicolor
でのfilter(Species == 'versicolor')のSpeciesを変数で指定する。
# 変数hoge hoge = 'Species' df %>% filter(hoge == 'versicolor') %>% head() # A tibble: 0 x 5 # … with 5 variables: Sepal.Length <dbl>, Sepal.Width <dbl>, Petal.Length <dbl>, # Petal.Width <dbl>, Species <fct>
このとき、変数hogeが評価されずに'hoge'列 が指定される。そんなものはないのでデータ抽出がされない。
expr(df %>% filter(hoge == 'versicolor') %>% head()) # df %>% filter(hoge == "versicolor") %>% head()
dplyr, ggplot, formula引数の書き方はNSE(Non-Standard Evaluation)という、表現式（expression）を与えて処理させる方法を取っている。
dplyrでのNSE
dplyrの表現式は例えば、
data %>% mutate(y = x + 1)
というように、x + 1は列名を直接書いていて、表現式で表している。これをSE(Standard Evaluation)という、関数に表現式ではなく値を与える形の場合は
data %>% mutate_(y =' x + 1')
という書き方になる。なお、mutate_のようにdplyrのNSE関数の一部は_を入れることでSE関数になる(非推奨)。
NSE関数の表現式では、上記の例で言えばx + 1は表現式となっているのでxは変数(シンボル)xとして通常使うことはできない。
この表現式はtidyevalのQuosure(quote(表現を評価しないままに留める） + closure(評価される環境を覚えておく）の造語)という特徴が影響しているそうだ。
そのため、表現式内において与える値はQuosureとして書く必要がある。また、そのままでは変数を与えていても評価しないままに留めるため、Quosureで渡した上でquoteを外す処理を外すことでxを変数xとして評価させることができる。
そのことを踏まえ、前述の以下の処理に戻ると、
# 変数hoge hoge = as.name('Species') #hogeをSpeciesというシンボルにする df %>% filter(hoge == 'versicolor') %>% # hogeはQuosure head() expr(df %>% filter(hoge == 'versicolor') %>% # hogeはQuosure head() ) # => df %>% filter(hoge == "versicolor") %>% head()
このfilter(hoge == 'versicolor')のfilter関数はNSE関数で、hogeはQuosure、つまりただの文字列hogeのため変数hogeとは別扱いになりhogeのまま扱われる。このQuosure内のhogeを変数としたい場合、hogeのquoteを外す!!を用いてこの部分は表現式のなかの文字列ではなく、変数だということ明示し評価対象とするとhogeは変数hogeとして評価される。
hoge = as.name('Species') #hogeをSpeciesというシンボルにする df %>% filter((!! hoge) == 'versicolor') %>% ## hogeというquoteとして渡したあと、!!でquoteを一時的に外して評価される head() # # A tibble: 6 x 5 # Sepal.Length Sepal.Width Petal.Length Petal.Width Species # <dbl> <dbl> <dbl> <dbl> <fct> # 1 7 3.2 4.7 1.4 versicolor # 2 6.4 3.2 4.5 1.5 versicolor # 3 6.9 3.1 4.9 1.5 versicolor # 4 5.5 2.3 4 1.3 versicolor # 5 6.5 2.8 4.6 1.5 versicolor # 6 5.7 2.8 4.5 1.3 versicolor expr(df %>% filter(!!hoge == 'versicolor') %>% # hogeはQuosure head() ) # =>df %>% filter(Species == "versicolor") %>% head()
ちなみに、hogeをシンボルにしない場合は評価はされるがfilter((!! hoge) == 'versicolor')はfilter(Species == 'versicolor')ではなくfilter('Species' == 'versicolor')になるのでデータ抽出はできない。
ggplotとformulaでのNSE
ggplotではaes()内に与えた値は、data$列名ではなく、与えたdataの列名を直接指定するような表現式になっている。
#NSE data %>% ggplot(aes(x = Sepal.Length, y = Sepal.Width)) # SE data %>% ggplot(aes(x = .$Sepal.Length, y = .$Sepal.Width)) #動かない
formulaでも同様に、dataの列名を直接指定する表現式としてNSEになっている。
# NSE data %>% lm(formula = Petal.Length ~ Petal.Width) # SE data %>% lm(formula = .$Petal.Length ~ .$Petal.Width) #動かない
自作関数の引数をNSEで渡す
自作関数の引数はシンボル(SE形式)なので一度NSE、つまり表現式に直さないといけない。そのため、表現式部分で変数を使いたい場合は一度enquo()でquosureに変換しないといけない。
df = as_tibble(iris) fun_select = function(df,x){ x = enquo(x) df %>% select(!!x) #df$Speciesの表現式になる } fun_select(df, Species)
fun_plot = function(df, x1, y1) { x1 = enquo(x1) y1 = enquo(y1) df %>% ggplot(aes(x = (!!x1), y =(!!y1))) + geom_point() } fun_plot(df, Sepal.Length, Sepal.Width)
Quosureはtidyevalという仕組みによって成り立っているみたいだが、NSEで書くと対応できるということ以外あまりそのあたりわかってないので詰まったらまた何か書くかも。
ちなみに、以下の過去記事で表面的な解決をしたが、根本的にはこの仕組みが原因っぽい。
参考
Programming with dplyr • dplyr
Rの関数定義でNSEを使う（表現式を引数にとれるようにする） | marketechlabo
↓このあたりのセルフ引用ツイート群 https://twitter.com/fronori/status/1002797245928955904
冒頭の疑問に戻って、Rで関数を作るときにexpressionを変数として認識（評価）させるtidyなやり方（tidyeval）は表現と環境をセットで表せるquosuresを用いる。 Quineのquasiquotationの概念が紹介され、実装される。まさに、意味深長。 https://t.co/jGJAtnJDO8 さっきの.dotsもpronounの一種だろう。
— Tetsuo Ishikawa (@fronori) 2018年6月2日
機械学習による予測確率は真の確率とは異なる
以下の記事では傾向スコアをロジスティク回帰で求めてその傾向スコアをもとにATEなどを求めた。
機械学習で出した確率は、予測確率が0.5未満ならラベル0、0.5以上ならラベル1にする、といったような分類器として使う場合は（おおむね）問題ないがそのとき出る予測確率そのものは真の確率とは異なるみたい。
そのため、機械学習を用いて確率自体を出したい場合は予測確率に対してProbability calibrationと呼ばれる補正が必要となる。
1.16. Probability calibration — scikit-learn 0.21.2 documentation
参考
傾向スコアによるマッチングを試す
傾向スコアによるマッチングを試す。内容は岩波DS3を、コードは以下を参考。
統計的因果推論(2): 傾向スコア(Propensity Score)の初歩をRで実践してみる - 渋谷駅前で働くデータサイエンティストのブログ
- 発売日: 2016/06/10
- メディア: 単行本（ソフトカバー）
傾向スコアとは
ランダム割当（RCT）となっていないような場合、「処置群になるかならないか」というのは各パネルの特徴（共変量）によって決まる傾向にある。この「施策群になる確率」を傾向スコアと呼ぶ。例えば、CM接触効果を測る場合、TVをよく見ている人ほどCM接触をしやすい（処置群になりやすい）ことから「TVをどれだけ見ているか」ということが共変量のひとつとなる。
そのため、処置群ダミーを目的変数とし、共変量を説明変数としてロジスティク回帰やプロビット回帰などをおこなうことで「処置群になるかどうかの確率（の推定値）」である傾向スコアを求めることができる。
この傾向スコアはあくまで傾向、さらにいえば推定値なので実際に施策がされたかどうかを示しているわけではないが、この傾向スコアによって補正をおこなうことで擬似的にRCT状態を組み込んだモデルを作成することができる。
なお、傾向スコアは共変量によって大きく変わるので、傾向スコアが妥当かどうかはc統計量（ROC曲線の下部面積）で判断でき、0.8以上であれば概ね問題ないとされているらしい。
また、確率を求めるためにロジスティック回帰以外も使うことができるが、機械学習によって求めた場合の確率はやや問題がありキャリブレーション（補正）が必要になるらしい。ここはそのうちちゃんと書きます。
統計的因果推論(4): 機械学習分類器による傾向スコアを調整してみる - 渋谷駅前で働くデータサイエンティストのブログ
ATE, ATT, ATU
母集団に対する施策の定量的な平均効果はATE(Average Treatment Effect:平均処置効果/因果効果)と呼ばれる。
また、処置群における平均処置効果はATT(Average Treatment effenct of the Treated)、対照群における平均処置効果はATU(Average Treatment effect of the Untreated)と呼ぶ。
ATTは「処置群に、施策を行った場合と行わなかった場合の差」となるため、その施策の効果測定をおこなうことができ、ATUは「対象群に、施策を行った場合と行わなかった場合の差」となるので、対照群に対して施策をおこなうかどうかの判断に使うことができる。
ATE, ATT, ATUの求め方
ATTは、「処置群に、施策を行った場合と行わなかった場合の差」から求めることができるが、施策を行った場合は行わなかった場合(あるいはその逆)の結果は得られない。そのため、「強い意味での無視可能性(strong ignorability)」という仮定をおくことで推定をおこなう。
この仮定は、
つまり「結果変数Yは共変量Xを条件付けると、介入群・対照群の割り付けZとは無関係（独立）で共変量Xのみに依存する」という仮定。この仮定を置くと、
となる。
は観測値から推定できないが、
は推定可能となるので、観測値だけから
を求めることができるため、
を求めることができる。
この仮定は例えば「Xのみに基づいて処置群に入れるかを判別する(k)」場合は成り立つし、「Yが高い人ほど処置群になりやすい」ような場合、つまりX以外がzを決めている場合は成り立たない。
ATEを求めるための手法は、処置群と対照群でマッチングさせたペアのY差平均を求める方法の他、岩波DS3では4つほど紹介されている。そのうち「IPW(Inverse Probability Weighting:逆確率重み付け法)」と「DR(Doubly Robust:二重にロバストな推定法)」の2つを考える。
IPWは「各々のサンプルに傾向スコアの逆数（対照群の場合は1から傾向スコアを引いた値）で重み付けをかけてATEを補正したもの」で、DRは「各々のサンプルに対してまず共変量による回帰モデルを作って再当てはめによって目的変数の予測値を出し、これと元の目的変数とで傾向スコアの逆数で加重平均させてATEを補正したもの」とのこと。計算式は岩波DS3のp.75参照。
ATE
サンプルデータlalondeを用いる。
まずは傾向スコアを求める
library(Matching) data(lalonde) # model作成 model <- glm(treat ~ age + educ + black + hisp + married + nodegr,lalonde,family=binomial) # fitさせ傾向スコアpsを求める ps <- model$fitted.values summary(model) # all: # glm(formula = treat ~ age + educ + black + hisp + married + nodegr, # family = binomial, data = lalonde) # # Deviance Residuals: # Min 1Q Median 3Q Max # -1.3941 -0.9933 -0.9237 1.3086 1.6633 # # Coefficients: # Estimate Std. Error z value Pr(>|z|) # (Intercept) 1.053028 1.047384 1.005 0.31471 # age 0.005917 0.014267 0.415 0.67833 # educ -0.063960 0.071354 -0.896 0.37005 # black -0.254369 0.363974 -0.699 0.48464 # hisp -0.829159 0.504230 -1.644 0.10009 # married 0.234241 0.266182 0.880 0.37886 # nodegr -0.838552 0.309383 -2.710 0.00672 ** # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # # (Dispersion parameter for binomial family taken to be 1) # # Null deviance: 604.20 on 444 degrees of freedom # Residual deviance: 589.43 on 438 degrees of freedom # AIC: 603.43 # # Number of Fisher Scoring iterations: 4 head(ps) # 1 2 3 4 5 6 # 0.4279364 0.2572819 0.5519755 0.3580844 0.4118540 0.3809886
c統計量はrms::lrmで求まる。
library(rms) model_lrm <- lrm(treat ~ age + educ + black + hisp + married + nodegr,lalonde) model_lrm # Logistic Regression Model # # lrm(formula = treat ~ age + educ + black + hisp + married + nodegr, # data = lalonde) # # Model Likelihood Discrimination Rank Discrim. # Ratio Test Indexes Indexes # Obs 445 LR chi2 14.77 R2 0.044 C 0.603 # 0 260 d.f. 6 g 0.407 Dxy 0.206 # 1 185 Pr(> chi2) 0.0221 gr 1.502 gamma 0.207 # max |deriv| 3e-12 gp 0.097 tau-a 0.100 # Brier 0.235
c統計量は0.603で基準値0.8に割と満たしてないがいったんこのままやる。
matchingで傾向スコアマッチングで処置群割当をおこない、グループ平均の差で平均処置効果を測定する。
# 共変量Xにpsを指定しているので傾向スコアが近い人がマッチングされる（傾向スコアマッチング）だが、ageとかを使うとageなどの属性をもとにしたマッチングとなる # match_ps <- Match(Y = lalonde$re78, Tr=lalonde$treat, X =ps) summary(match_ps) # Estimate... 2302 # AI SE...... 764.97 # T-stat..... 3.0093 # p.val...... 0.0026187 # # Original number of observations.............. 445 # Original number of treated obs............... 185 # Matched number of observations............... 185 # Matched number of observations (unweighted). 492
傾向スコアマッチング推定だと、ATEは2302。これは対照群/処置群のマッチングペア毎にYの差分を出し、その平均値を出している。
次に、IPW推定を試す。
# IPW y <- lalonde$re78 z1 <- lalonde$treat ipwe1 <- sum((z1*y)/ps)/sum(z1/ps) ipwe0 <- sum(((1-z1)*y)/(1-ps))/sum((1-z1)/(1-ps)) ipwe1 - ipwe0 # 1632.463
ATEのIPW推定量は1632と傾向マッチング推定の2302より小さく出た。
また、直接周辺期待値を計算しないで、傾向スコアの逆数で重み付け（補正）をおこなって回帰してtreatの効果をみることでも計算は可能。
この場合は、周辺期待値だけでなく標準誤差なども合わせて求められる。
このときのweightは
式は以下のp.408
http://www.cardio.med.tohoku.ac.jp/newmember/pdf/ms/hmbr46-399.pdf
y <- lalonde$re78 ivec1 <- lalonde$treat # treatダミー ivec0 = rep(1, nrow(lalonde))-ivec1 # untreatダミー ivec <- cbind(ivec1,ivec0) iestp1 <- z1/ps # treatの場合の逆確率重み付け(untreatの場合0) iestp0 <- (1-z1)/(1-ps) # untreatの場合の逆確率重み付け(treatの場合0) iestp <- iestp1+iestp0 # 片方は0になるのでtreat,untreatを区別せず逆確率重み付け値 ipwe_gs <- lm(y ~ ivec + 0, weights=iestp) # 逆確率重み付けをした加重最小二乗法。ivec1(treatダミー)=0の場合は0*重み=0になる summary(ipwe_gs) # Call: # lm(formula = y ~ ivec + 0, weights = iestp) # # Weighted Residuals: # Min 1Q Median 3Q Max # -10534 -5886 -2389 4072 90256 # # Coefficients: # Estimate Std. Error t value Pr(>|t|) # ivecivec1 6196.0 454.7 13.63 <2e-16 *** # ivecivec0 4563.5 454.1 10.05 <2e-16 *** # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # # Residual standard error: 9584 on 443 degrees of freedom # Multiple R-squared: 0.3929, Adjusted R-squared: 0.3902 # F-statistic: 143.4 on 2 and 443 DF, p-value: < 2.2e-16
6196.0 - 4563.5 = 1632.5となり、直接求めた場合と同じ(丸めているので下2桁は不明だが）。
次にDR推定。
# DR dre <- function(data, target, treat, ps, formula) { n <- nrow(data) y <- data[target] data1 <- data[data[treat]==1,] data0 <- data[data[treat]==0,] model1 <- lm(formula=formula, data=data1) model0 <- lm(formula=formula, data=data0) fitted1 <- predict(model1, data) fitted0 <- predict(model0, data) dre1 <- (1/n)*sum(y+((z1-ps)/ps)*(y-fitted1)) dre0 <- (1/n)*sum(((1-z1)*y)/(1-ps)+(1-(1-z1)/(1-ps))*fitted0) return(c(dre1, dre0)) } ret <- dre(lalonde, "re78", "treat", ps, + re78 ~ age + educ + black + hisp + married + nodegr) ret # 6185.915 4566.401 ret[1] - ret[2] # 1619.514
ATEのDR推定量は1620、IPW推定量の1632、傾向マッチング推定の2302よりも小さい。
処置群の中には処置群なのに傾向スコアが低い人、逆に対照群なのに傾向スコアが高い人がいる。傾向スコアによる逆確率重み付けはこういう状況を補正する。そのため単に処置群/対照群の平均差を出した傾向マッチング推定と、傾向スコアで補正をおこなったDR推定量、IPW推定量で結果が大きく異なっている。
なお、CBPSパッケージに様々な関数があるとのこと。
ATT
引き続き、「処置群に、施策を行った場合と行わなかった場合の差」であるATTを求める。
このときのweightは
式はATE同様に以下のp.409
http://www.cardio.med.tohoku.ac.jp/newmember/pdf/ms/hmbr46-399.pdf
コードが微妙にわかるようなわからんような。。。
model <- glm(treat ~ age + educ + black + hisp + married + nodegr,lalonde,family=binomial) ps <- model$fitted.values # fitさせ傾向スコアpsを求める ivec1 <- lalonde$treat # treatダミー iestp1_ATT <- ivec1 #treatの場合の重み（そのまま） ivec0 <- rep(1, nrow(lalonde)) - lalonde$treat # untreatダミー iestp0_ATT <- ivec0 * ps / (1-ps) #untreatの場合の重み ivec <- cbind(ivec1,ivec0) iestp_ATT <- iestp1_ATT + iestp0_ATT ipwe_ATT_re74 = lm(lalonde$re74 ~ ivec + 0, weights=iestp_ATT) summary(ipwe_ATT_re74) # Call: # lm(formula = lalonde$re74 ~ ivec + 0, weights = iestp_ATT) # # Weighted Residuals: # Min 1Q Median 3Q Max # -2853 -2096 -1753 -1180 32945 # # Coefficients: # Estimate Std. Error t value Pr(>|t|) # ivecivec1 2095.6 363.4 5.766 1.52e-08 *** # ivecivec0 2226.0 362.9 6.134 1.90e-09 *** # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # # Residual standard error: 4943 on 443 degrees of freedom # Multiple R-squared: 0.1379, Adjusted R-squared: 0.134 # F-statistic: 35.44 on 2 and 443 DF, p-value: 5.29e-15
結果より、2095.6 - 2226.0 ≒ -130から処置することで130ほど負の効果がある。
・・・なんかおかしい気がする。傾向スコアのc統計量は0.603だったからか？
参考
調査観察データにおける因果推論(2) - 傾向スコアとIPW推定量，二重にロバストな推定量 - About connecting the dots.
調査観察データにおける因果推論(3) - Rによる傾向スコア，IPW推定量，二重にロバストな推定量の算出 - About connecting the dots.
統計的因果推論に関するスライドとRのサンプルコード | かものはしの分析ブログ
差分の差分法(DID)を試す
差分の差分法(DID)の勉強のために以下の記事を参考にする。
https://fisproject.jp/2016/05/difference-in-differences-using-r/
使用するデータは 、A-Gまでの7ヶ国についてy, y_bin, x1, x2, x3, opinion を1990年から10年間に渡り記録したパネルデータとのこと。
require(foreign) library(tidyverse) # Panel data panel = read.dta("https://dss.princeton.edu/training/Panel101.dta") head(panel) # country year y y_bin x1 x2 x3 opinion # 1 A 1990 1342787840 1 0.2779036 -1.1079559 0.28255358 Str agree # 2 A 1991 -1899660544 0 0.3206847 -0.9487200 0.49253848 Disag # 3 A 1992 -11234363 0 0.3634657 -0.7894840 0.70252335 Disag # 4 A 1993 2645775360 1 0.2461440 -0.8855330 -0.09439092 Disag # 5 A 1994 3008334848 1 0.4246230 -0.7297683 0.94613063 Disag # 6 A 1995 3229574144 1 0.4772141 -0.7232460 1.02968037 Str agree panel %>% ggplot(aes(year, y, color = country)) + geom_line() + geom_vline(xintercept = 1994 ,colour="black")
1994年にE, F, G国で施策が実施されたとして, この3ヶ国を処置群, それ以外を対照群として前後の y の変化をみるため, ダミー変数を追加する。
このとき追加するダミー変数は、
- 処置群ダミー(treat)
- 施策前後ダミー(postperiod)
- DID効果ダミー(処置群ダミー * 施策前後ダミー = 施策の効果がある期間ダミー)
の3つを入れ、これらを説明変数とし、以下のOLSを行えばよい。
なお、Xは処置群・施策前後以外でYに影響を与える要因、eは誤差を表す。
この重回帰で何故DID効果が求められるかというと、簡易化のため で考えると、
- Treatで、処置群/対照群の差（グループ差）でのYに与えている影響のうち時点(PostPre)に依らない部分を統制
- PostPreで、時点のY差に与えている影響のうちグループ差(PostPre)に依らない部分を統制
- Treat * PostPreで、グループ差・時点差の両方によるYの差 = 処置群で施策後の効果 = DID効果
となっている。
統計的因果推論(1): 差分の差分法(Difference-in-Differences)をRで回してみる - 渋谷駅前で働くデータサイエンティストのブログ
こちらのサイトの図
でいえば、A-C(E-D)の差がTreatの係数、A-E(C-D)の差がPostPreの係数、B-EがTreat * PostPreの係数で表される(※A-Cなどの-は引き算という意味ではない)。
B = C(切片) + (A-C Treat差) + (A-E PostPre差) + (B-E DID差) となっていることからもわかる。
なお、Xを置かない場合の仮定として処置群と対照群で施策が無い場合同じような動きをするという「平行トレンド」の仮定が必要となる(Xを置く場合はXで調整も可能)。
今回のデータは施策実施前の1994年以前をみるとまぁ概ねどれも同じような傾向の動きをしているのでとりあえずOKとしておく。
この仮定がおけるクラスタの作り方がわからない場合、傾向スコアマッチングなどをおこなうようだがこれは以降の記事で試す。
# 処置群ダミー panel$treated = ifelse(panel$country == "E" | panel$country == "F" | panel$country == "G", 1, 0) # 施策前後ダミー panel$postperiod = ifelse(panel$year >= 1994, 1, 0) # DID効果ダミー panel$did = panel$postperiod * panel$treated # difference in differences model.did = lm(y ~ treated + postperiod + did, data = panel) # => # p-value 0.0882 : The effect is significant at 10% with the treatment having a negative effect. summary(model.did) # Call: # lm(formula = y ~ treated + time + did, data = panel) # # Residuals: # Min 1Q Median 3Q Max # -9.768e+09 -1.623e+09 1.167e+08 1.393e+09 6.807e+09 # # Coefficients: # Estimate Std. Error t value Pr(>|t|) # (Intercept) 3.581e+08 7.382e+08 0.485 0.6292 # treated 1.776e+09 1.128e+09 1.575 0.1200 # time 2.289e+09 9.530e+08 2.402 0.0191 * # did -2.520e+09 1.456e+09 -1.731 0.0882 . # --- # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 # # Residual standard error: 2.953e+09 on 66 degrees of freedom # Multiple R-squared: 0.08273, Adjusted R-squared: 0.04104 # F-statistic: 1.984 on 3 and 66 DF, p-value: 0.1249
結果から、DIDの効果は10%水準では-2.520e+09有意に負の影響がある（駄目じゃん）となっている。
参考
差分の差分分析（Difference-in-differences design） – 医療政策学×医療経済学
- 作者: 山本勲
- 出版社/メーカー: 中央経済社
- 発売日: 2015/10/21
- メディア: 単行本
- この商品を含むブログ (1件) を見る
各モデルで使われる変数重要度についてまとめる
はじめに
機械学習タスクにおいて、ブラックボックス性を減らし解釈性を上げるためにどの特徴量がどれくらい効いたかを示す必要性が増えている。また、特徴量選択やフィーチャーエンジニアリングの際にも参考になる。そのため、機械学習タスクをやるにあたって必須のスキルとなっているのでまとめてみる。
なお、構成や大まかな流れは以下のスライドからかなり引用しています(というか、下記スライドに個人的なメモを追加して書いただけに近い)。 speakerdeck.com
機械学習モデルに対しての特徴量解釈は主に以下の3つの観点がある、
- どの特徴量が重要か
- 各特徴量が予測にどう影響するか
- ある予測結果に対して特徴量がどう寄与するか
各観点によって使う手法が変わり、以下でそれぞれの内容をまとめる。
どの特徴量が重要か
モデルが重要視している要因がわかる。
Feature Importance
モデルに使用する特徴量の重要度を数値化したもの。
モデルに合わせたアルゴリズム(Model Dependent)と、モデルに関わらないアルゴリズム（Model Independent）がある。
実際の使い方としては、意思決定においてどの特徴量に対して優先して改善を図るか、といった際に使用したり、Feature Enginiaringにおいてとりあえず特徴量を全部ぶち込んでから各特徴量を使うかどうかの指標にしたりする。
Model Dependent
- Random Forest:評価指標をより改善する特徴量が重要
- 勾配Boosting系(XGBoost , LightGBM, CatBoost...)
- Gain, cover, weight(frequency)
…など。
どの変数を分割したら上記判別基準がどう変わったかで重要度を計算する。 ちなみに、Scikit-learnのRandom ForestはGini imuprityベース。 それぞれの詳細は以下。
Random Forest
Random Forestで計算できる特徴量の重要度 - なにメモ
ランダムフォレスト — 苦労する遊び人の玩具箱 1 ドキュメント
www.slideshare.net
ランダムフォレスト 特徴量の重要度（C++の実装例つき） - じじいのプログラミング
XGboost
Xgboost : feature_importanceのimportance_type算出方法 - Qiita
Gradient Boosted Tree (Xgboost) の取り扱い説明書 - Qiita
Model Independent アルゴリズムがモデルに依存しないので作成した各モデルに対して統一的に重要度を見ることができる。
- Permutation Important:ある特徴量をランダムに交換したとき、予測精度が大きく落ちる場合、重要度が高いと判定される。
Permutation Importance | Kaggle
Python: 特徴量の重要度を Permutation Importance で計測する - CUBE SUGAR CONTAINER
permutation importance 〜特徴量の重要度の測り方〜 - TJ Media Studies
https://bunseki-train.com/permutation_importance_and_eli5_1/
変数重要度とPartial Dependence Plotでブラックボックスモデルを解釈する - Dropout
各特徴量が予測にどのように影響しているか
ある特徴量を変化させたときに予測がどう変化するかという関係性がわかる。
例えば、特徴量X1の値が3を境にYへの貢献度が一気に増える、など。
実際の使い方としては、Feature Importanceを見たあとに、重要な特徴量のいくつかに対して改善を図ろうと考える。その際に例えば特徴量X1は値が変わってもあまり貢献度は変わらないが特徴量X2は値を変えると大きく貢献度が変わるという場合、特徴量X2に焦点を当てて改善していこう、といったことが考えられる。
Partial Dependence
特徴量を変化させたときの出力の平均的な変化。 Partial Dependence Plot(PDP)で可視化できる。 ただし、特徴量同士の相関が強い場合は信用できない。
平均ではなく、各レコードについて個別に関係を見ていくIndividual Conditional Expectation Plot(ICE plot)というものもある。
PDPでは、変数間の交互作用の影響がわからないがICE plotでは確認ができる。
5.1 Partial Dependence Plot (PDP) | Interpretable Machine Learning
Partial Dependence Plots • pdp
5.2 Individual Conditional Expectation (ICE) | Interpretable Machine Learning
機械学習モデルを解釈する方法 Permutation Importance / Partial Dependence Plot - 子供の落書き帳 Renaissance
ある予測結果に対して特徴量がどう寄与するか
何故その予測値だったかわかる。例えば、予測値Y1では特徴量X1が大きく寄与していたけど、予測値Y2では特徴量X2の寄与が大きかった、など。
SHAP, LIME
説明したい予測値に対してシンプルな局所的な分類器をつくり、以下の条件から貢献度を考える
- 1.局所的正確性:データ点の直上では値が一致
- 2.値がゼロの変数の無影響性:データ点で値が0の変数は貢献度0
- 3.複数モデル比較での無矛盾性:簡単化の前後でモデル間の貢献度の大小関係が同じ
SHAP(SHapley Additive exPlanation) value
条件を全て満たす貢献度の組。
- お互いの値次第で出力が変わるモデル
- 最終的な出力への貢献を公平に変数に分配
複数の予測値の傾向の可視化も可能
Feature ImportanceやPartial Dependenceに近い使い方も可能
5.9 Shapley Values | Interpretable Machine Learning
機械学習モデルを解釈する指標SHAPについて – 戦略コンサルで働くデータサイエンティストのブログ
続・機械学習モデルを解釈する方法 SHAP value - 子供の落書き帳 Renaissance
LIME(Local Interpretable Model-agnostic Explanation)
最適化問題で貢献度を算出
5.7 Local Surrogate (LIME) | Interpretable Machine Learning
機械学習モデルの予測結果を説明するための力が欲しいか...？ - クソして寝ろ
LIMEでXGBoostとLightGBMの結果を解釈する - Qiita
Grad-CAM
画像系に使う？画像系は使う機会がなさそうなのであまりちゃんと調べてないです。
Class Active Mapping(CAM) Grad-CAM:勾配ベースのCAM DNNに対して、各判定クラスで予測への影響が強い部分を可視化
SHAP valueを使ってみる
SHAP valueを実際に試してみる。コードは以下の記事のものを拝借。
なお、データはkaggleのHouse Pricingコンペデータをテキトーに加工して作っている。SalePrice(対数化済)を目的変数として、簡易化のため変数重要度がTOP5を試す。
また、モデルはRandomForestを使う。
import numpy as np # linear algebra import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv) # for visualization import matplotlib.pyplot as plt import seaborn as sns # import data train = pd.read_csv(‘../input/train.csv’) test = pd.read_csv(‘../input/test.csv’) # string label to categorical values # 全objectタイプをlabel encode from sklearn.preprocessing import LabelEncoder for i in range(train.shape[1]): if train.iloc[:,i].dtypes == object: lbl = LabelEncoder() lbl.fit(list(train.iloc[:,i].values) + list(test.iloc[:,i].values)) #ここにある項目で変換する train.iloc[:,i] = lbl.transform(list(train.iloc[:,i].values)) test.iloc[:,i] = lbl.transform(list(test.iloc[:,i].values)) # keep ID for submission train_ID = train[‘Id’] test_ID = test[‘Id’] # split data for training y_train = train[‘SalePrice’] X_train = train.drop([‘Id’,‘SalePrice’], axis=1) X_test = test.drop(‘Id’, axis=1) # dealing with missing data #missing多い場合はdrop、少ない場合はmean X_train = X_train.drop([‘LotFrontage’,‘MasVnrArea’,‘GarageYrBlt’], axis=1) X_train = X_train.fillna(X_train.median()) X_test = X_test.drop([‘LotFrontage’,‘MasVnrArea’,‘GarageYrBlt’], axis=1) X_test = X_test.fillna(X_train.median()) # add SF X_train[‘TotalSF’] = X_train[‘TotalBsmtSF’]+X_train[‘1stFlrSF’]+X_train[‘2ndFlrSF’] X_test[‘TotalSF’] = X_test[‘TotalBsmtSF’]+X_test[‘1stFlrSF’]+X_test[‘2ndFlrSF’] y_train = np.log(y_train)
scikit-learnのRandomForest.feature_importances_
まずはじめに、一般的によくみる、giniベースのアルゴリズムであるscikit-learnのfeature_importances_関数を用いて変数重要度を見る。
# feature importance using random forest from sklearn.ensemble import RandomForestRegressor rf = RandomForestRegressor(n_estimators=80, max_features=‘auto’) rf.fit(X_train, y_train) print(‘Training done using Random Forest’) ranking = np.argsort(-rf.feature_importances_) #降順で配列のindexを返す # use the top 5 features only X_train = X_train.iloc[:,ranking[:5]] X_test = X_test.iloc[:,ranking[:5]]
SHAP
SHAPで各特徴量が予測にどう効いたかを見てみる。
まずは、各特徴量の値に対しての、予測への影響度であるSHAP値を算出。
import shap shap.initjs() # jupyterの表示用 explainer = shap.TreeExplainer(rf) shap_values = explainer.shap_values(X_train)
SHAP 要約(平均)
summary_plotで各特徴量の寄与度絶対値平均を図示。
これはginiベースのものと概ね同じになった。 ただし、「なにを重要とするか」という指標が各アルゴリズムで異なるため、必ずしも同じような値になるわけではない。
shap.summary_plot(shap_values, X_train, plot_type="bar")
SHAP 要約(個々)
summary_plotで各特徴量の、個々の予測値(点1つ）に対する寄与度を図示。色は特徴量の値の大きさ。
shap.summary_plot(shap_values, X_train)
インスタンス（各レコード）が点となっていて、値の大小で色づけされている。OverrallQualをみてみると貢献度(SHAP value)が上がるほど色が大きく（赤）なっていることからキレイに相関が出ているし、他でも概ねその傾向がある。
ある予測値に対しての貢献度の内訳
force_plotで、あるインスタンス(以下のコードでは1番目のレコード)に対する各特徴量の予測値に対する寄与の内訳を表している。赤が予測値を増加させる特徴量、青が減少させる特徴量となる。base valueで示されている数値は予測値の期待値。黒太字の数値は予測値。棒の長さはどの程度の寄与をしたかの絶対値、となる。
shap.force_plot(explainer.expected_value, shap_values[0,:], X_train.iloc[0,:])
実際の値で例示すると、赤の合計は大体0.23程度で青の合計が0.04程度の長さとなっているので、12.21(output value) ≒ 12.02(base value) + 0.23(赤) - 0.04(青)になっていることを考えればわかりやすい。
個々で見ていくと、OverrallQualの長さは0.17程度だが、これはOverrallQualが7になることでSalePrice(対数化済)が0.17程度増加した、と考えることができる。
全予測値に対してのある特徴量と出力値の関係
先程はforce_plotである1つのインスタンスについての内訳を見たが今回は全インスタンスを一気に見る。
以下の図ではYearBuiltとoutput value(予測値)について出力している。図にカーソルを合わせると個々の内訳が見れる。また、x,y軸は任意のものに変えることができる。
shap.force_plot(base_value=explainer.expected_value, shap_values=shap_values, features=X_train)
色の赤/青は、YearBuiltがある値のときにYearBuiltが正に寄与したか負に寄与したかを表している。同じYearBuilt値なのに赤/青が混在している部分では正に寄与したり負に寄与したりとなっているのは、他の特徴量の影響によるもの。
全予測値に対してのある特徴量の貢献度の内訳
ある特徴量の値とSHAP valueの関係。
shap.dependence_plot(ind="YearBuilt",interaction_index='YearBuilt', shap_values=shap_values, features=X_train)
上図をみると、値が大きくなるほどSHAP valueも増加していることがわかる。
また、特徴量同士の関係もみることができる。
shap.dependence_plot(ind='YearBuilt',interaction_index='TotalSF', shap_values=shap_values, features=X_train)
上手では縦軸にYearBuiltのSHAP value、横軸にYearBuilt、色としてTotalSFの値が表示されている。YearBuiltの値が変わってもTotalSFの大小の傾向はあまり変わらない。ただ、なんとなく1950頃を境に、1950以前はTotalSFが高いとYearBuiltのSHAP valueは高くなっているが、以降ではotalSFが低いとYearBuiltのSHAP valueは高くなっているようにも見える。
参考
https://orizuru.io/blog/machine-learning/shap/
複数のモデルを管理する
purrrとbroomの使い方をマスターするために。以下の記事をトレースする。 内容自体は引用元記事の方がちゃんとしているのでそちらを読んでもらいたい。本記事は、読むにあたって理解が薄い人（=自分)用の補足メモを書きながらのトレースとなる。そのため、本筋外の解説が多く、ごちゃごちゃしていることに留意してもらいたい。
purrrとbroomで複数の回帰モデルを効率的に管理する - Dropout
引用記事にもあるが、複数のモデルを作成する際に個々でオブジェクトを作っていくとデータなどに変更があった際に全てに忘れずに適用しなければいけなくなったり、どのオブジェクトがどの係数になっているかという係数用オブジェクトを各モデル毎に更にオブジェクトで作っていくと更にややこしいことになる。また、モデルの比較も面倒。
そのため、各モデルを1つの計算処理でおこなえるようにしつつ、結果を1つのオブジェクトで一元管理ができるならば非常に便利。
これを実現するために、purrr::map()とbroom::tidy(), broom::glance()を用いる。
library(tidyverse) library(tidymodels) library(patchwork) # 可視化用 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()) df = diamonds
データの構造を見る
目的変数を正規分布にする
今回、予測モデルとしてダイヤモンドのpriceを、carat, clarity(透明度ランク)を用いてOLSで推定する。
その際、（本論のモデル管理とは別だが）、OLSを用いるためpriceの確率分布を（できる限り）正規分布にするためにpriceに対数を取って推定する必要がある。
# 正規分布にする df %>% ggplot(aes(log(price))) + geom_histogram(fill = "#56B4E9", color = "white") + theme_minimal2
log(price)とcaratは非線形関係なので、caratにも対数をとり線形にする必要がある。
# 線形回帰のために線形関係にする # カラットvs値段 # カラットが上がるほど値段があがる df %>% sample_frac(0.1) %>% #データが多いので減らす ggplot(aes(log(carat), log(price))) + geom_point(color = "#0072B2", alpha = 0.5) + theme_scatter
直感的には透明度が上がると値段が上がりそうだが、データをみると逆の関係になっている。
これは、透明度が上がるほどカラット数が小さくなる...おそらく、カラット数が大きく透明度が高いダイヤは作りづらい（もしくは需要が少ないからあまり作らない）ということが反映されている。
そのため、カラット数が交絡となり、値段と透明度に負の相関が現れる結果となっている。
# 透明度vs値段 # => 透明度が上がるほど値段がさがる g1 = df %>% ggplot(aes(clarity, log(price))) + geom_boxplot(fill = "#56B4E9") + theme_minimal2 # 透明度vsカラット # =>透明度が上がるほどカラットが下がる g2 = df %>% ggplot(aes(clarity, log(carat))) + geom_boxplot(fill = "#56B4E9") + theme_minimal2 g1 | g2
モデリング
使用するモデルの作成
上記の結果、
- 価格、透明度は対数を取る
- 交絡を考慮し、モデルにはカラット数を入れる
このことを踏まえ、検証できるように以下の3モデルを作成する。
- 1.価格の対数を、透明度のみで説明するモデル
- 2.価格の対数を、透明度とカラット数で説明するモデル
- 3.価格の対数を、透明度とカラット数の対数で説明するモデル
# 回帰用前処理 df_input = df %>% mutate_if(is.ordered, factor, ordered = FALSE)# factor(ordered = FALSE)に変換 formulas = c(log(price) ~ clarity, log(price) ~ clarity + carat, log(price) ~ clarity + log(carat)) %>% enframe("model_no", "formula") #vactor to DF # model_no formula # <int> <list> # 1 1 <S3: formula> # 2 2 <S3: formula> # 3 3 <S3: formula>
enframe()
enframe(x, name = "name", value = "value")はvector/listをtibbleに変換する関数。上記コードでは、パイプでxを渡し、引数1,2はそれぞれmodel_no, formulaを指定している。これは列名に対応している（デフォルトだと、name, valueになる）。
以下はhelpのサンプルコード
enframe(1:3) # => # # A tibble: 3 x 2 # name value # <int> <int> # 1 1 1 # 2 2 2 # 3 3 3 enframe(c(a = 5, b = 7)) # => # # A tibble: 2 x 2 # name value # <chr> <dbl> # 1 a 5 # 2 b 7 enframe(list(one = 1, two = 2:3, three = 4:6)) # => # # A tibble: 3 x 2 # name value # <chr> <list> # 1 one <dbl [1]> # 2 two <int [2]> # 3 three <int [3]>
回帰モデルの結果の整理
先程作成したformulasをmap()を用いてデータにfitさせる。
# モデルとその評価を管理しやすい形にする df_result = formulas %>% mutate(model = map(formula, lm, data = df_input), # lm(data = df_input, 各formulas$formula) tidied = map(model, tidy), # 回帰モデルの係数 to tidy tibble glanced = map(model, glance)) # 回帰モデルの評価 # => # # 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]>
map()
map(.x, .f, ...)関数を用いると、xのそれぞれのデータに対して関数/formula fを適応させる。また、f以降の引数はfの引数として扱われる。
今回の処理をmap関数を用いないで、formulasの1行目(log(price) ~ clarity)だけに適応させると以下になる。
lm(data = df_input, log(price) ~ clarity) # => # Call: # lm(formula = log(price) ~ clarity, data = df_input) # # Coefficients: # (Intercept) claritySI2 claritySI1 clarityVS2 clarityVS1 clarityVVS2 clarityVVS1 # 8.0276 0.1392 -0.1797 -0.2646 -0.3029 -0.4968 -0.7047 # clarityIF # -0.6225
そのため、model列はformulasの各行に入っているformulaに対して適応させた結果が入っている。
tidy()
tidy(x, ...)はxのオブジェクトをtidyな形のtibbleに変換する。今回の場合は、model列の内容(係数)をxとして変換している。
今回のtidied列に入る結果もmap関数を用いない場合以下のようになる。
lm(data = df_input, log(price) ~ clarity) %>% #model列 tidy() # modelの結果の係数をtidyなtibble # A tibble: 8 x 5 # term estimate std.error statistic p.value # <chr> <dbl> <dbl> <dbl> <dbl> # 1 (Intercept) 8.03 0.0363 221. 0. # 2 claritySI2 0.139 0.0377 3.69 2.28e- 4 # 3 claritySI1 -0.180 0.0373 -4.81 1.49e- 6 # 4 clarityVS2 -0.265 0.0374 -7.08 1.50e-12 # 5 clarityVS1 -0.303 0.0379 -7.99 1.41e-15 # 6 clarityVVS2 -0.497 0.0389 -12.8 2.44e-37 # 7 clarityVVS1 -0.705 0.0398 -17.7 7.20e-70 # 8 clarityIF -0.622 0.0432 -14.4 5.03e-47
glance
glance(x, ...)はxのmodelオブジェクトのsummaryを1行に変換する。今回の場合は、model列の内容をxとして変換している。
今回のglanced列に入る結果もmap関数を用いない場合以下のようになる。
lm(data = df_input, log(price) ~ clarity) %>% #model列 glance() # modelの結果のモデル評価(summaryの一部)を1行で # => # # A tibble: 1 x 11 # r.squared adj.r.squared sigma statistic p.value df logLik AIC BIC deviance # <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> # 1 0.0511 0.0510 0.988 415. 0 8 -75907. 1.52e5 1.52e5 52694. # # … with 1 more variable: df.residual <int>
モデルの係数比較
先程の結果df_resultの係数部分(tidied列)の中身を見るために`unnest()'する。
# model結果の係数を下二桁 df_coef = df_result %>% select(model_no, tidied) %>% unnest() %>% mutate_if(is.double, round, digits=2) # => # # 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
このままでは比較がしづらいので、横持ちに変換する。
# 係数を横持ち化 df_coef %>% mutate(term = fct_inorder(term)) %>% # デフォルト（アルファベット順）から出てきた順にする dplyr::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
結果を見ると、
- mode1では、カラット数の交絡の影響で、透明度が上がるほど大きな負の係数となる
- model2では、カラット数の交絡の影響を統制したため、透明度が上がるほど大きな正の係数となる。
- model3では、model2とほぼ同様だがcaratに対数を取って統制をとったことで透明度の影響（係数）にやや変化が出ている。
modelの評価
モデルの評価が入っているglanced列に対しても同様に比較する。
# モデルの評価 df_result %>% dplyr::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 AIC BIC # <int> <dbl> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> # 1 1 0.05 0.05 0.99 415. 0 8 -75907. 1.52e5 1.52e5 # 2 2 0.87 0.87 0.37 43737. 0 9 -23024. 4.61e4 4.62e4 # 3 3 0.97 0.97 0.19 187918. 0 9 13378. -2.67e4 -2.66e4 # # … with 2 more variables: deviance <dbl>, df.residual <int>
このとき、自由度調整済み決定係数adj.r.squaredを見るとmodel3が1番性能が高くなっている。model2と比較すると、OLSを用いているので価格とカラット数を線形関係に変換した方が性能が高くなる、というあ妥当な結果になっている。
まとめ
今回のコードをまとめると以下のようになる。
流れとしては、 - 1.回帰モデルを1オブジェクト(formulas)に作成する - 2.map関数を用いることで、formulasに対しfit, 結果, 評価を整理して1オブジェクトに格納 - 3.見たい部分に対してunnest()で結果を見る
# モデル比較の下準備------- # 回帰用前処理 df_input = df %>% mutate_if(is.ordered, factor, ordered = FALSE)# factor(ordered = FALSE)に変換 # 回帰モデル # 1.価格を透明度のみで説明するモデル # 2.透明度とカラット数で説明するモデル # 3.上と同じだが、カラット数に対数をとったモデル formulas = c(log(price) ~ clarity, log(price) ~ clarity + carat, log(price) ~ clarity + log(carat)) %>% enframe("model_no", "formula") #vactor to DF # モデルとその評価を管理しやすい形にする df_result = formulas %>% mutate(model = map(formula, lm, data = df_input), # lm(data = df_input, 各formulas$formula) tidied = map(model, tidy), # 回帰モデルの係数 to tidy tibble glanced = map(model, glance)) # 回帰モデルの評価 # モデルの係数比較---------- # model結果の係数を下二桁 df_coef = df_result %>% dplyr::select(model_no, tidied) %>% unnest() %>% mutate_if(is.double, round, digits=2) # 係数を横持ち化 df_coef %>% mutate(term = fct_inorder(term)) %>% # デフォルト（アルファベット順）から出てきた順にする dplyr::select(model_no, term, estimate) %>% spread(model_no, estimate) # モデルの評価比較------- # model結果の評価を下二桁 df_result %>% dplyr::select(model_no, glanced) %>% unnest() %>% mutate_if(is.double, round, digits=2)
