{mlogit}パッケージ(多項ロジットモデル)の学習記録

 

【学習動機】

少し複雑なロジットモデルをRで推定するために, {mlogit}パッケージの使い方を覚えたい.
調べていたら, Train先生が{mlogit}の使い方のエクササイズを書いているのを発見した.
Kenneth Train’s exercises using the mlogit package for R
http://www.idg.pl/mirrors/CRAN/web/packages/mlogit/vignettes/Exercises.pdf
これに沿って学習していく. 今回は, 多項ロジットモデルの節を扱う.
ドキュメントの全てを扱うわけではない. また変数名やコードなど自分の理解しやすい形に変えているところもあるので注意されたい.  

【追記:2019/06/16】  

Exercise 1: Multinomial logit model  

ググって探索していたら, 二週間前にこれがアップロードされていたことに気付いた. 英語が読めるならこっちを読んだ方が絶対にいい気がする.  

 

【学習内容】

1. 多項ロジットモデル

library(mlogit)
library(tidyverse)
ライブラリ読み込み 
 

カリフォルニア州の900家計の暖房システムの選択データ
新築でエアコンの設置された戸建て住宅
選択肢は5つ:gc (gas central), gr (gas room), ec (electric central), er (electric room), hp (heat pump)

data("Heating",package = "mlogit")
head(Heating)
##   idcase depvar  ic.gc  ic.gr  ic.ec  ic.er   ic.hp  oc.gc  oc.gr  oc.ec
## 1      1     gc 866.00 962.64 859.90 995.76 1135.50 199.69 151.72 553.34
## 2      2     gc 727.93 758.89 796.82 894.69  968.90 168.66 168.66 520.24
## 3      3     gc 599.48 783.05 719.86 900.11 1048.30 165.58 137.80 439.06
## 4      4     er 835.17 793.06 761.25 831.04 1048.70 180.88 147.14 483.00
## 5      5     er 755.59 846.29 858.86 985.64  883.05 174.91 138.90 404.41
## 6      6     gc 666.11 841.71 693.74 862.56  859.18 135.67 140.97 398.22
##    oc.er  oc.hp income agehed rooms region
## 1 505.60 237.88      7     25     6 ncostl
## 2 486.49 199.19      5     60     5 scostl
## 3 404.74 171.47      4     65     2 ncostl
## 4 425.22 222.95      2     50     4 scostl
## 5 389.52 178.49      2     25     6 valley
## 6 371.04 209.27      6     65     7 scostl

変数は以下のようになっている.
idcase:id
depvar:暖房システム
ic.alt:暖房システムaltの設置コスト
oc.alt:暖房システムaltの年間運用コスト
pb.alt:比率oc.alt / ic.alt
income:世帯の年収
agehed:世帯主の年齢
rooms:部屋の数
region:ncostl(北部沿岸地域)、scostl(南部沿岸地域)、mountn(山岳地帯)、valley(中央渓谷地域)

選択の割合を確認

table(Heating$depvar)/900    #900人で割って割合を出す
## 
##         gc         gr         ec         er         hp 
## 0.63666667 0.14333333 0.07111111 0.09333333 0.05555556

gc : 64%, gr : 14%, ec : 7%, er : 9%, hp : 6%

wide形式のデータになっているので, mlogit.data()で{mlogit}パッケージで推定できるように形式を変更.

H_long <- mlogit.data(Heating, shape = "wide", choice = "depvar",
                       varying = c(3:12))

各家計の選択関係なく平均

H_long %>% group_by(alt) %>% summarise(ave_income = mean(income),
                                       ave_age = mean(agehed),
                                       ave_rooms = mean(rooms),
                                       ave_ic = mean(ic),
                                       ave_oc = mean(oc))
## # A tibble: 5 x 6
##   alt   ave_income ave_age ave_rooms ave_ic ave_oc
##   <chr>      <dbl>   <dbl>     <dbl>  <dbl>  <dbl>
## 1 ec          4.64    42.9      4.42   825.   477.
## 2 er          4.64    42.9      4.42   984.   430.
## 3 gc          4.64    42.9      4.42   777.   172.
## 4 gr          4.64    42.9      4.42   922.   154.
## 5 hp          4.64    42.9      4.42  1046.   219.

選択を条件付けて平均

H_long %>% filter(depvar == TRUE) %>% 
  group_by(alt) %>% summarise(ave_income = mean(income),
                              ave_age = mean(agehed),
                              ave_rooms = mean(rooms),
                              ave_ic = mean(ic),
                              ave_oc = mean(oc))
## # A tibble: 5 x 6
##   alt   ave_income ave_age ave_rooms ave_ic ave_oc
##   <chr>      <dbl>   <dbl>     <dbl>  <dbl>  <dbl>
## 1 ec          4.70    44.3      4.59   789.   451.
## 2 er          4.62    39.3      4.52   928.   407.
## 3 gc          4.68    43.3      4.40   781.   173.
## 4 gr          4.36    43.9      4.37   933.   155.
## 5 hp          4.88    40.5      4.46  1027.   213.

ベースモデル

設置コストと年間運用コストで説明するモデル. コストなので効用を下げるはずなので, 係数は負になると予想される.

model_0 <- mlogit(depvar ~ ic + oc | 0, data = H_long)
summary(model_0)
## 
## Call:
## mlogit(formula = depvar ~ ic + oc | 0, data = H_long, method = "nr")
## 
## Frequencies of alternatives:
##       ec       er       gc       gr       hp 
## 0.071111 0.093333 0.636667 0.143333 0.055556 
## 
## nr method
## 4 iterations, 0h:0m:0s 
## g'(-H)^-1g = 1.56E-07 
## gradient close to zero 
## 
## Coefficients :
##       Estimate  Std. Error z-value  Pr(>|z|)    
## ic -0.00623187  0.00035277 -17.665 < 2.2e-16 ***
## oc -0.00458008  0.00032216 -14.217 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -1095.2

全選択肢の確率

model_0 %>% fitted(outcome=FALSE) %>% head()
##           ec         er        gc        gr         hp
## 1 0.09545811 0.05094155 0.4642482 0.3166757 0.07267644
## 2 0.05831236 0.03698376 0.4482688 0.3696122 0.08682283
## 3 0.08396029 0.03195210 0.6220979 0.2250557 0.03693402
## 4 0.11457049 0.09663155 0.2883760 0.4375632 0.06285874
## 5 0.07692237 0.03737185 0.4188408 0.2806727 0.18619232
## 6 0.13353597 0.05281490 0.5279799 0.1725121 0.11315710

900人の平均選択確率

model_0 %>% fitted(outcome=FALSE) %>% apply(2, mean)
##         ec         er         gc         gr         hp 
## 0.10413057 0.05141477 0.51695653 0.24030898 0.08718915

データの割合と比較

table(Heating$depvar)/900
## 
##         gc         gr         ec         er         hp 
## 0.63666667 0.14333333 0.07111111 0.09333333 0.05555556

フィット的に見て, あまりいいモデルとは言えなそう.

coef(model_0)["oc"]/coef(model_0)["ic"]
##        oc 
## 0.7349453

年間運用コストを1ドル削減するのに支払える設置コストを0.73ドル多く支払えるという結果になっている.
年間運用コストは毎年払うものであるのに対して, 設置コストは1回のみ. この結果は合理的でない.
やはりモデルが良くなさそう.

生涯コストモデル

割引率を使って年間運用コストと設置コストを生涯コストに変換したモデル.

#H_long <- H_long %>% mutate(lcc = ic + oc/0.12)
#ここでdplyr::mutateを使うとmlogit.data型が消えてしまうので使えない
H_long$lcc <- H_long$ic + H_long$oc/0.12
model_lcc <- mlogit(depvar ~ lcc | 0, data = H_long)
summary(model_lcc)
## 
## Call:
## mlogit(formula = depvar ~ lcc | 0, data = H_long, method = "nr")
## 
## Frequencies of alternatives:
##       ec       er       gc       gr       hp 
## 0.071111 0.093333 0.636667 0.143333 0.055556 
## 
## nr method
## 5 iterations, 0h:0m:0s 
## g'(-H)^-1g = 9.32E-05 
## successive function values within tolerance limits 
## 
## Coefficients :
##        Estimate  Std. Error z-value  Pr(>|z|)    
## lcc -7.1585e-04  4.2761e-05 -16.741 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -1248.7

ベースモデルとの比較で, 尤度比検定を行う.

lrtest(model_0,model_lcc)
## Likelihood ratio test
## 
## Model 1: depvar ~ ic + oc | 0
## Model 2: depvar ~ lcc | 0
##   #Df  LogLik Df  Chisq Pr(>Chisq)    
## 1   2 -1095.2                         
## 2   1 -1248.7 -1 306.93  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#年間運用コストの係数と設置コストの係数との関係が制約条件となっている
qchisq(0.05, df = 1, lower.tail = FALSE)
## [1] 3.841459

model_0が支持される.

選択肢別切片モデル

各選択ごとに異なる切片を設定したモデル. hpの選択肢を0として基準にしている.

model_con <- mlogit(depvar ~ ic + oc, data = H_long, reflevel = "hp")
summary(model_con)
## 
## Call:
## mlogit(formula = depvar ~ ic + oc, data = H_long, reflevel = "hp", 
##     method = "nr")
## 
## Frequencies of alternatives:
##       hp       ec       er       gc       gr 
## 0.055556 0.071111 0.093333 0.636667 0.143333 
## 
## nr method
## 6 iterations, 0h:0m:0s 
## g'(-H)^-1g = 9.58E-06 
## successive function values within tolerance limits 
## 
## Coefficients :
##                   Estimate  Std. Error z-value  Pr(>|z|)    
## ec:(intercept)  1.65884594  0.44841936  3.6993 0.0002162 ***
## er:(intercept)  1.85343697  0.36195509  5.1206 3.045e-07 ***
## gc:(intercept)  1.71097930  0.22674214  7.5459 4.485e-14 ***
## gr:(intercept)  0.30826328  0.20659222  1.4921 0.1356640    
## ic             -0.00153315  0.00062086 -2.4694 0.0135333 *  
## oc             -0.00699637  0.00155408 -4.5019 6.734e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -1008.2
## McFadden R^2:  0.013691 
## Likelihood ratio test : chisq = 27.99 (p.value = 8.3572e-07)

このモデルから得られる900人の平均選択確率

model_con %>% fitted(outcome=FALSE) %>% apply(2, mean)
##         hp         ec         er         gc         gr 
## 0.05555556 0.07111111 0.09333333 0.63666667 0.14333333

データの割合と比較

table(Heating$depvar)/900
## 
##         gc         gr         ec         er         hp 
## 0.63666667 0.14333333 0.07111111 0.09333333 0.05555556

データでの割合とモデルで予測される平均選択確率が一致していてフィットは良さそう.

支払い意思額

coef(model_con)["oc"]/coef(model_con)["ic"]
##       oc 
## 4.563385

割引率

1/(coef(model_con)["oc"]/coef(model_con)["ic"])
##        oc 
## 0.2191356

割引率が22%と計算でき, 前のモデルに比べたら妥当.

収入あたり設置コストモデル

設置コストを収入あたりにしたモデル.

model_income <- mlogit(depvar ~ oc + I(ic/income), data = H_long)
summary(model_income)
## 
## Call:
## mlogit(formula = depvar ~ oc + I(ic/income), data = H_long, method = "nr")
## 
## Frequencies of alternatives:
##       ec       er       gc       gr       hp 
## 0.071111 0.093333 0.636667 0.143333 0.055556 
## 
## nr method
## 6 iterations, 0h:0m:0s 
## g'(-H)^-1g = 1.03E-05 
## successive function values within tolerance limits 
## 
## Coefficients :
##                  Estimate Std. Error z-value  Pr(>|z|)    
## er:(intercept)  0.0639934  0.1944893  0.3290  0.742131    
## gc:(intercept)  0.0563481  0.4650251  0.1212  0.903555    
## gr:(intercept) -1.4653063  0.5033845 -2.9109  0.003604 ** 
## hp:(intercept) -1.8700773  0.4364248 -4.2850 1.827e-05 ***
## oc             -0.0071066  0.0015518 -4.5797 4.657e-06 ***
## I(ic/income)   -0.0027658  0.0018944 -1.4600  0.144298    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -1010.2
## McFadden R^2:  0.011765 
## Likelihood ratio test : chisq = 24.052 (p.value = 5.9854e-06)

対数尤度が選択肢別切片モデルより低くなっている.

 

model_income2 <- mlogit(depvar ~ oc + ic | income, data = H_long, reflevel = "hp")
summary(model_income2)
## 
## Call:
## mlogit(formula = depvar ~ oc + ic | income, data = H_long, reflevel = "hp", 
##     method = "nr")
## 
## Frequencies of alternatives:
##       hp       ec       er       gc       gr 
## 0.055556 0.071111 0.093333 0.636667 0.143333 
## 
## nr method
## 6 iterations, 0h:0m:0s 
## g'(-H)^-1g = 6.27E-06 
## successive function values within tolerance limits 
## 
## Coefficients :
##                   Estimate  Std. Error z-value  Pr(>|z|)    
## ec:(intercept)  1.95445797  0.70353833  2.7780 0.0054688 ** 
## er:(intercept)  2.30560852  0.62390478  3.6954 0.0002195 ***
## gc:(intercept)  2.05517018  0.48639682  4.2253 2.386e-05 ***
## gr:(intercept)  1.14158139  0.51828845  2.2026 0.0276231 *  
## oc             -0.00696000  0.00155383 -4.4792 7.491e-06 ***
## ic             -0.00153534  0.00062251 -2.4664 0.0136486 *  
## ec:income      -0.06362917  0.11329865 -0.5616 0.5743846    
## er:income      -0.09685787  0.10755423 -0.9005 0.3678281    
## gc:income      -0.07178917  0.08878777 -0.8085 0.4187752    
## gr:income      -0.17981159  0.10012691 -1.7958 0.0725205 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Log-Likelihood: -1005.9
## McFadden R^2:  0.01598 
## Likelihood ratio test : chisq = 32.67 (p.value = 1.2134e-05)
#尤度比検定
lrtest(model_con, model_income2)
## Likelihood ratio test
## 
## Model 1: depvar ~ ic + oc
## Model 2: depvar ~ oc + ic | income
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1   6 -1008.2                     
## 2  10 -1005.9  4 4.6803     0.3217
#ワルド検定
waldtest(model_con, model_income2)
## Wald test
## 
## Model 1: depvar ~ ic + oc
## Model 2: depvar ~ oc + ic | income
##   Res.Df Df  Chisq Pr(>Chisq)
## 1    894                     
## 2    890  4 4.6456     0.3256
#スコア検定
scoretest(model_con, model_income2)
## 
##  score test
## 
## data:  depvar ~ oc + ic | income
## chisq = 4.6761, df = 4, p-value = 0.3222
## alternative hypothesis: unconstrained model

収入は選択に影響しないのかもしれない.

ヒートポンプの設置コスト割引するケース

hp(ヒートポンプ)の設置費用の選択に対して払い戻しを10%行うとき, シェアはどのようになるか予測する.

X <- model.matrix(model_con)
head(X)
##      ec:(intercept) er:(intercept) gc:(intercept) gr:(intercept)      ic
## 1.ec              1              0              0              0  859.90
## 1.er              0              1              0              0  995.76
## 1.gc              0              0              1              0  866.00
## 1.gr              0              0              0              1  962.64
## 1.hp              0              0              0              0 1135.50
## 2.ec              1              0              0              0  796.82
##          oc
## 1.ec 553.34
## 1.er 505.60
## 1.gc 199.69
## 1.gr 151.72
## 1.hp 237.88
## 2.ec 520.24
alt <- index(H_long)$alt
chid <- index(H_long)$chid
#hpの設置コストを1割引に
X[alt == "hp", "ic"] <- X[alt == "hp", "ic"] * 0.9
eXb <- as.numeric(exp(X %*% coef(model_con)))
SeXb <- tapply(eXb, chid, sum)
P <- eXb/SeXb[chid]
P <- matrix(P, ncol = 5, byrow = TRUE)
P_pred <- apply(P, 2, mean) %>% t()
colnames(P_pred) <- c("ec","er","gc","gr","hp")
P_pred 
##              ec         er        gc        gr        hp
## [1,] 0.07045486 0.09247026 0.6306444 0.1419681 0.0644623

割引前と比較

model_con %>% fitted(outcome=FALSE) %>% apply(2, mean)
##         hp         ec         er         gc         gr 
## 0.05555556 0.07111111 0.09333333 0.63666667 0.14333333

割引によって, ヒートポンプの割合が5.5%から6.4%に上昇することがこのモデルから予測される.

【学習予定】

次は, ネスティッドロジットモデルの節を扱う.

ロジットモデルの人工データ生成とOVBの学習記録

 

【学習動機】

傾向スコアをそのうち仕事で使いそうな雰囲気があるので、傾向スコアを学び直しておこう、というかちゃんと学ぼうってことになったわけなんだが、その前にロジットモデルを学び直しておこうと思った。そんなこんなで、Trainの“Discrete Choice Methods with Simulation”をパラパラと読んでいる。

https://www.amazon.co.jp/Discrete-Choice-Methods-Simulation-Kenneth/dp/0521747384

ロジットモデルを学んですぐのチュートリアルでは、実データでやることが多いが、それらは本当に推定できているのかという疑問を持った。
今回はシミュレーションデータを作りながら、ちゃんと推定できるのか確認していき、学び直す。ついでに、ロジットモデルでもOVB(Omitted Variable Bias : 欠落変数バイアス)が生じることをシミュレーションで確認する。

 

【学習内容】

ロジットモデルに関して、ざっくりとメモ。離散選択モデルはランダム効用モデルから導出される。個人は選択によって効用が最大になるような選択肢を選ぶものとし、効用(U)を代表的効用(V)と誤差項を分解し、誤差項にガンベル分布を当てはめると、選択肢iの選択確率がexp(Vi)/sum(exp(Vi))となってロジットモデルが導かれる。Vを線形モデルで当てはめ、選択確率の積に対数をとって対数尤度関数を導出し、最尤法でもって、パラメータを推定する。対数尤度関数は線形のパラメータに関して大局的に凹であることはMcFadden(1974)で証明されているようだ(標準的なロジットモデルではなく、ネスティッドロジットモデルになるとそうではないとTrainのテキストのGEVの章に書かれている)。

以下でロジットモデルの人工データを生成して、ロジットモデルで推していく。

ライブラリ読み込み。

library(ggplot2)    #visualization

関数を先に定義しておく。

#ロジット関数
logit <- function(p){
  return( log(p) - log(1-p) )
}
#ロジスティク関数
logistic <- function(y){
  return( 1 / (1 + exp(-y)))
}

モデルの説明

真のモデルのパラメータ設定\(\beta_{1}=-4, \beta_{2}=8\)でサンプルサイズは1000, シミュレーション回数は1000としておく。

# true parameter
beta1 <- -4
beta2 <- 8
# sample size
n <- 1000
# # of simulation
m <- 1000
#seed
set.seed(1234)

説明変数x1とx2を作っておく。図示もしておく。

x1 <- c(runif(n = n/2, min = 0, max = 800), runif(n = n/2, min = 400, max = 1600))/1000
x2 <- rbinom(n = n, size = 1, prob = 0.6)
plot(x1)

plot(x2)

今回の真のモデルについての説明。
選択確率は以下のように、
\(\beta_{1}=-4, \beta_{2}=8\)がパラメータとなり、
\(y^{*} = x1*\beta_{1} + x2*\beta_{2}\)
\(P = \frac{1}{1+exp(-y^{*})}\)とする。
選択は以下のように、
\(y\) ~ \(Bernoulli(P)\)
データが発生するものとする。

このデータに対して、\(y^{*} = x1*\beta_{1} + x2*\beta_{2}\)で真のモデルに等しいモデル(model0と呼ぶ)で推定する場合と、\(y^{*} = x1*\beta_{1}\)でx2が脱落している過少定式のモデル(modelmisと呼ぶ)で推定する場合を見ていく。

err <- evd::rgumbel(n = n, loc = 0, scale = 1)
ystar <- x1*beta1 + x2*beta2
utility <- ystar + err
p <- logistic(ystar)
y <- rbinom(n = n, size = 1, prob = p)
  
df <- data.frame(
    ystar = ystar,
    utility = utility,
    p = p,
    y = y,
    x1 = x1,
    x2 = x2
  )
  
model0 <- glm(formula = y ~ 0 + x1 + x2, family = binomial(link = "logit"), data = df)
modelmis <- glm(formula = y ~ 0 + x1, family = binomial(link = "logit"), data = df)

summary(model0)
## 
## Call:
## glm(formula = y ~ 0 + x1 + x2, family = binomial(link = "logit"), 
##     data = df)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -3.11148  -0.18597   0.06014   0.12358   2.68804  
## 
## Coefficients:
##    Estimate Std. Error z value Pr(>|z|)    
## x1  -3.8564     0.3471  -11.11   <2e-16 ***
## x2   7.8470     0.5388   14.56   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1386.29  on 1000  degrees of freedom
## Residual deviance:  323.78  on  998  degrees of freedom
## AIC: 327.78
## 
## Number of Fisher Scoring iterations: 7
summary(modelmis)
## 
## Call:
## glm(formula = y ~ 0 + x1, family = binomial(link = "logit"), 
##     data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.5064  -1.2777   0.9916   1.0831   1.1771  
## 
## Coefficients:
##    Estimate Std. Error z value Pr(>|z|)    
## x1  0.46689    0.07947   5.875 4.23e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1386.3  on 1000  degrees of freedom
## Residual deviance: 1350.5  on  999  degrees of freedom
## AIC: 1352.5
## 
## Number of Fisher Scoring iterations: 4

当然、model0はよく推定できているし、modelmisよりAICが低い。modelmisの\(\beta_{1}\)の推定値が全然違うし、符号も逆という結果になっている。

ちなみに、今回は使わないが、{evd}パッケージのrgumbel()関数でガンベル分布からの乱数を発生させられる。
utility出してるけど、ここでは、なんの意味もない。utilityのmax取る選択肢をとる実験を何回も繰り返すと、選択確率がざっくり出てくるだろうけど。

シミュレーション

上の人工データ生成と推定を1000回繰り返す。model0で推定されたパラメータはb1_0とb2_0に保存し、modelmisで推定されたパラメータはb1_misに保存する。

#for estimates
b1_0 <- numeric(m)
b1_mis <- numeric(m)
b2_0 <- numeric(m)

#seed
set.seed(1234)

for(i in 1:m){
  err <- evd::rgumbel(n = n, loc = 0, scale = 1)
  ystar <- x1*beta1 + x2*beta2
  utility <- ystar + err
  p <- logistic(ystar)
  y <- rbinom(n = n, size = 1, prob = p)
  
  df <- data.frame(
    ystar = ystar,
    utility = utility,
    p = p,
    y = y,
    x1 = x1,
    x2 = x2
  )
  
  model0 <- glm(formula = y ~ 0 + x1 + x2, family = binomial(link = "logit"), data = df)
  modelmis <- glm(formula = y ~ 0 + x1, family = binomial(link = "logit"), data = df)
  
  b1_0[i] <- model0$coefficients[[1]]
  b1_mis[i] <- modelmis$coefficients[[1]]
  b2_0[i] <- model0$coefficients[[2]]
}
df_est <- data.frame(
  b1_0 = b1_0,
  b1_mis = b1_mis,
  b2_0 = b2_0
)

#histogram of estimates of beta1
p0 <- ggplot(df_est,aes(x=b1_0)) + geom_histogram(colour="grey") 
pmis <- ggplot(df_est,aes(x=b1_mis)) + geom_histogram(colour="grey") 

#b1
gridExtra::grid.arrange(p0, pmis)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#b2
ggplot(df_est,aes(x=b2_0)) + geom_histogram(colour="grey") 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

model0は真のモデルと定式が同じだからしっかり推定できている。一方、modelmisは先ほど確認した通り、\(\beta_{1}\)の符号から違うという結果に。

失敗例

推定が上手くいった例を上に載せたが, 失敗例も役立つかもしれないので載せておく. 失敗例を2つ. 映画のエンディングのNGシーン集的な.

サンプルサイズが100のケース

nを減らすと正しく収束しなかった. サンプルサイズを100として同じシミュレーションをやる

# sample size
n <- 100
set.seed(1234)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: glm.fit: fitted probabilities numerically 0 or 1 occurredと警告が出る.
推定が上手くいってない. 失敗というか上手くいってない推定が何回か起きてる感じか.

x1を1000で割らないケース

サンプルサイズを元の戻して, x1を1000で割らずに同じシミュレーションをする.

# sample size
n <- 1000
set.seed(1234)
x1 <- c(runif(n = n/2, min = 0, max = 800), runif(n = n/2, min = 400, max = 1600))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

最適化の問題か. 探索範囲とかそういう問題で推定が上手くいってないと思われる. 最適化の実践的な問題に詳しくないので, 事実かどうかはよく分からない.

【学習予定】

次はネスティッドロジットモデルあたりのシミュレーションデータ作りながら確認か, 傾向スコアの方をやってみたい.
ロジットモデルのOVBの理論が簡単に読めるならそれも読んでおきたい.

{easyPubMed}パッケージの学習記録

 

【学習動機】

PubMedをRで扱えないものかと思ってググってみたら, {easyPubMed}というパッケージがあった. 検索のときに, ある程度キーワードで絞った後, Rでabstractとかの単語検索して自分の探したい文献を見つけ出したり, 情報を洗い出すのに役立つだろうなと思ったので, {easyPubMed}の使い方を学んでいく.

過去の日本語の記事で{easyPubMed}を扱ったものが一件あったが, パッケージのアップデートによって関数や色々変わっていて, そのままでは使えない状態になっている.

www.karada-good.net

そのため, 公式ドキュメントを眺めながら学習していくことにして, その記録を残す.  

https://cran.r-project.org/web/packages/easyPubMed/easyPubMed.pdf

【学習内容】

私はPubMedについて全く詳しくないのでPubMedの説明もここではしない. PubMedの使い方の基本は東大の医学部図書館などが公開している.  

https://www.lib.m.u-tokyo.ac.jp/manual/pubmedmanual.pdf

以下, {easyPubMed}パッケージの使い方を, ドキュメントを見ながら確認していく.

今回の検索の事例として, 最近, 尿路結石が気になるので尿路結石について調べていく. 尿路結石の治療として外から衝撃波で壊すという話を周りの尿路結石経験者が語っていたのを覚えていた. 体外衝撃波結石破砕術というものらしい. ガイドラインによると,  

「30 数年前の内視鏡を用いた体内式結石破砕術(TULとPNL)と,その数年後に導入された体外衝撃波結石破砕術(ESWL)の開発により,腎切石術,腎盂切石術,尿管切石術などの開放手術が激減し,コペルニクス的な転回を遂げたことは広く知られるところである。」とある. そんなわけで, ESWLのメタ解析をしてそうな文献を探し出すことを今回の事例とする.

尿路結石症診療ガイドライン 2013年版 | Mindsガイドラインライブラリ

ライブラリ読み込み.

library(easyPubMed)
library(tidyverse)    #宇宙
library(stringr)    #文字列の処理
library(lubridate)    #時間の計算

尿路結石をMeSHで, とりあえず見てみる. 検索数を見る. get_pubmed_ids()はリストで簡易的な情報が返ってくる.

#尿路結石("Urinary Calculi")
Query1 <- '"urinary calculi"[MeSH]'
pubmedidlist <- get_pubmed_ids(Query1)
pubmedidlist$Count    #ヒット数
## [1] "34458"

体外衝撃波結石破砕術(ESWL : Extracorporeal Shock Wave Lithotripsy)とメタ解析で絞っていく.

Query2 <- '"extracorporeal shock wave lithotripsy"'
Query3 <- 'meta-analysis'
Query <- str_c('(',Query1,Query2,Query3,')')
Query
## [1] "(\"urinary calculi\"[MeSH]\"extracorporeal shock wave lithotripsy\"meta-analysis)"
get_pubmed_ids(Query)$Count
## [1] "28"

28まで絞れたが, 新しめの知見を得たいので10年以内の文献に絞る.

researchdate <- "20190609"    #検索実行日の決定
tenyearsago <- ymd(researchdate) - years(10) + days(3)
tenyearsago <- tenyearsago %>% as.character() %>% str_replace_all("-","/")
researchdate <- ymd(researchdate) %>% as.character() %>% str_replace_all("-","/")
Querytime <- str_c('AND("',tenyearsago,'"[PDat] : "',researchdate, '"[PDat])' )
Query <- Query %>% str_c(Querytime)
pubmedidlist <- get_pubmed_ids(Query)
pubmedidlist$Count
## [1] "20"

day(3)とずらしているのは, PubMedで10年以内のフィルタをかけると, クエリは’AND (“2009/06/12”[PDat] : “2019/06/09”[PDat])’となる. 10年の期間の計算は何をしているのかは正しく理解しているわけではないが, このへん5年だと12日が11日になるところを見ると, おそらく, 365*10(or 5)を単純に引いて+1日しているもので, うるう年のうるう日分だけ前に来ると考えられるので, 以上のような計算を行った.
結果, 20まで絞れた.
次に, fetch_all_pubmed_idf()で, その20個のpubmedidを手に入れて, カンマでつなげてurlを作成すると, 該当pubmedidのurlができる.

pubmedid <- fetch_all_pubmed_ids(pubmedidlist)
pubmedid
##  [1] "29686172" "28068364" "27798945" "26530230" "26383613" "26382788"
##  [7] "26352272" "26211003" "25909212" "25860144" "25681251" "25449204"
## [13] "24919112" "24863324" "24162890" "22592707" "22161396" "21855945"
## [19] "19889063" "19560860"
pubmedidstr <- pubmedid %>% str_c(collapse = ",")
pubmedidstr %>% str_c("https://www.ncbi.nlm.nih.gov/pubmed/", .)
## [1] "https://www.ncbi.nlm.nih.gov/pubmed/29686172,28068364,27798945,26530230,26383613,26382788,26352272,26211003,25909212,25860144,25681251,25449204,24919112,24863324,24162890,22592707,22161396,21855945,19889063,19560860"

www.ncbi.nlm.nih.gov

fetch_pubmed_dataで詳しい情報入手. XMLで返ってくる.

papers <- fetch_pubmed_data(pubmedidlist, encoding = "ASCII")
papers %>% custom_grep("ArticleTitle", "char") %>% length()
## [1] 20
papers %>% custom_grep("Abstract", "char") %>% length()
## [1] 19

custom_grep()はabstractがない文献(本やコメントなど)がNAで入らないので, 上のように件数が上手く合わなくて個人的に使いにくい.
というわけで, ごり押し気味だけど以下の方法で情報を抽出する.
articles_to_list()で各文献の情報がリストで手に入る. そのリストの中身をそれぞれ, article_to_df()でデータフレームで情報(PMID, DOI, article title, article abstract, publication date (year, month, day), journal name (title, abbreviation), keywords, and a set of author-specific info (names, affiliation, email address))が手に入る. max_charsはabstract全文は10000あれば取れるだろうという設定. 著者が複数いると別のオブザベーションとして入るので, ここでは, 欲しい情報の抜き出しには1行目だけ使う.

articlelist <- papers %>% articles_to_list()
narticles <- length(articlelist)
pmid <- numeric(narticles)
year <- numeric(narticles)
title <- character(narticles)
abstract <- character(narticles)

for(i in 1:narticles){
  df <- articlelist[[i]] %>% article_to_df(max_chars = 10000)
  pmid[i] <- df[1,]$pmid %>% as.numeric()
  year[i] <- df[1,]$year %>% as.numeric()
  title[i] <- df[1,]$title
  abstract[i] <- df[1,]$abstract
}

result <- data.frame(
  pmid = pmid,
  year = year,
  title = title,
  abstract = abstract
)

これでidと出版年とタイトルとアブストラクトの情報をもったデータフレームが手に入った.

アブストラクトに, 体内式結石破砕術(TUL, PNL)のどちらかが入ってるの探す.

result %>% mutate(flg_TUL = str_detect(abstract,"TUL"),
                  flg_PNL = str_detect(abstract, "PNL")) %>% 
  filter(flg_TUL==TRUE | flg_PNL) %>% select(-abstract)
##       pmid year
## 1 25449204 2015
##                                                                                                                                                                                       title
## 1 Systematic review and meta-analysis of the clinical effectiveness of shock wave lithotripsy, retrograde intrarenal surgery, and percutaneous nephrolithotomy for lower-pole renal stones.
##   flg_TUL flg_PNL
## 1   FALSE    TRUE

という感じで, それらしいものが一つだけ得られた.

www.ncbi.nlm.nih.gov

【学習予定】

PubMedの使い方としてはお粗末だったかもしれないが, {easyPubMed}に関してはこんな感じで使える. 今後は未定. 

RでNaNやNAやInfをデータフレーム全体で置換するときの覚書

 

【用途】
RでデータをいじっていてInfやNaN, NAを何か置換するとき, 気付かないうちに, 文字だと思っていた変数が数値に変わっていることがある. そんなときのための覚書.

【内容】
文字で入れていたはずの変数が中で因子型(factor)になっていると, 置換の操作をしたときに因子を表す数字に変換されていることがあったので, 型をしっかり把握しておくこと. 因子型(factor)から文字型(charactor)にしておけば, 変換されずに済む.

ライブラリ読み込み.

library(tidyverse)
## Warning: package 'tidyverse' was built under R version 3.5.3
## -- Attaching packages ------------------------------------------------------------------------------- tidyverse 1.2.1 --
## √ ggplot2 3.0.0     √ purrr   0.2.5
## √ tibble  1.4.2     √ dplyr   0.7.6
## √ tidyr   0.8.1     √ stringr 1.3.1
## √ readr   1.1.1     √ forcats 0.3.0
## -- Conflicts ---------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

具体例のデータフレーム. tibbleにした方が型が見えやすいのでtibbleにしておく.

df <- data.frame(a = c("焼肉", "寿司", "二郎"),
                 b = c(Inf,1,0),
                 c = c(3,3,NaN),
                 d = c(2,NA,1)) %>% as_tibble()

df
## # A tibble: 3 x 4
##   a         b     c     d
##   <fct> <dbl> <dbl> <dbl>
## 1 焼肉    Inf     3     2
## 2 寿司      1     3    NA
## 3 二郎      0   NaN     1

このデータフレームのInfやNaN,NAを0に置き換えたいとする.

df %>% mutate_all(funs(ifelse(is.infinite(.),0,.))) %>% 
  mutate_all(funs(ifelse(is.nan(.),0,.))) %>% 
  mutate_all(funs(ifelse(is.na(.),0,.)))
## # A tibble: 3 x 4
##       a     b     c     d
##   <int> <dbl> <dbl> <dbl>
## 1     2     0     3     2
## 2     1     1     3     0
## 3     3     0     0     1

とやってみると, 焼肉も寿司も二郎もなくなってしまった.

一個ずつ変換を試みる. さらに追加で2種類変換の実験.

df %>% mutate_all(funs(ifelse(is.infinite(.),0,.)))
## # A tibble: 3 x 4
##       a     b     c     d
##   <int> <dbl> <dbl> <dbl>
## 1     2     0     3     2
## 2     1     1     3    NA
## 3     3     0   NaN     1
df %>% mutate_all(funs(ifelse(is.nan(.),0,.)))
## # A tibble: 3 x 4
##       a     b     c     d
##   <int> <dbl> <dbl> <dbl>
## 1     2   Inf     3     2
## 2     1     1     3    NA
## 3     3     0     0     1
df %>% mutate_all(funs(ifelse(is.na(.),0,.)))
## # A tibble: 3 x 4
##       a     b     c     d
##   <int> <dbl> <dbl> <dbl>
## 1     2   Inf     3     2
## 2     1     1     3     0
## 3     3     0     0     1
df %>% replace_na(list(b=0,c=0,d=0))
## # A tibble: 3 x 4
##   a         b     c     d
##   <fct> <dbl> <dbl> <dbl>
## 1 焼肉    Inf     3     2
## 2 寿司      1     3     0
## 3 二郎      0     0     1
df %>% na_if(0)
## # A tibble: 3 x 4
##   a         b     c     d
##   <fct> <dbl> <dbl> <dbl>
## 1 焼肉    Inf     3     2
## 2 寿司      1     3    NA
## 3 二郎     NA   NaN     1

最初の3つは焼肉寿司二郎が2,1,3に置き換わってしまっている.
naに関してはreplace_naが使えるが, リストによる指定が面倒くさい.

因子型となっているのが悪そうなので, 文字型に変換してやってみる.

df$a <- as.character(df$a)

df %>% mutate_all(funs(ifelse(is.infinite(.),0,.)))
## # A tibble: 3 x 4
##   a         b     c     d
##   <chr> <dbl> <dbl> <dbl>
## 1 焼肉      0     3     2
## 2 寿司      1     3    NA
## 3 二郎      0   NaN     1
df %>% mutate_all(funs(ifelse(is.nan(.),0,.)))
## # A tibble: 3 x 4
##   a         b     c     d
##   <chr> <dbl> <dbl> <dbl>
## 1 焼肉    Inf     3     2
## 2 寿司      1     3    NA
## 3 二郎      0     0     1
df %>% mutate_all(funs(ifelse(is.na(.),0,.)))
## # A tibble: 3 x 4
##   a         b     c     d
##   <chr> <dbl> <dbl> <dbl>
## 1 焼肉    Inf     3     2
## 2 寿司      1     3     0
## 3 二郎      0     0     1
df %>% replace_na(list(b=0,c=0,d=0))
## # A tibble: 3 x 4
##   a         b     c     d
##   <chr> <dbl> <dbl> <dbl>
## 1 焼肉    Inf     3     2
## 2 寿司      1     3     0
## 3 二郎      0     0     1
df %>% na_if(0)
## # A tibble: 3 x 4
##   a         b     c     d
##   <chr> <dbl> <dbl> <dbl>
## 1 焼肉    Inf     3     2
## 2 寿司      1     3    NA
## 3 二郎     NA   NaN     1

上手くいっているので全部やると

df %>% mutate_all(funs(ifelse(is.infinite(.),0,.))) %>% 
  mutate_all(funs(ifelse(is.nan(.),0,.))) %>% 
  mutate_all(funs(ifelse(is.na(.),0,.)))
## # A tibble: 3 x 4
##   a         b     c     d
##   <chr> <dbl> <dbl> <dbl>
## 1 焼肉      0     3     2
## 2 寿司      1     3     0
## 3 二郎      0     0     1

なんとか上手くいった. どういうメカニズムでこの問題が起きているか分かってないうえ, もっといい置換方法があるよう気もするが, 場当たり的にとりあえずこれで解決したので書き留めておく.

【記憶の検索キーワード】
NA, NaN, Inf, 置換, 全部, 文字, 数値, 因子, R

Sports Analyst Meetup #2の後日譚々

Sports Analyst Meetup #2@データスタジアムにてLTを行ってきました。

そのスライドがこちらになります。

 

speakerdeck.com

 

自分の発表力不足もあり、きちんと言いたいことが伝わったかどうかは分かりませんが、もし詳細が気になった方がいらっしゃいましたら、以下の記事を読んでみたり、お気軽にお声をかけてください。

コードや詳細はこちらの記事にまとめてあります。

moratoriamuo.hatenablog.com

この分析は、遺伝的アルゴリズムだとかガウス過程とか高度なものや面白いデータは一切使っておらず、言わばサイコロを振るだけの簡単なシミュレーションなのですが、簡易なモデル設定から本質を突けている面白いものだと、自分では思っています。

モデルの発想元は大学に入ってから頭を使ってバドミントンをするようになったところにあります。中学高校とバドミントンをしていましたが、当時はどんな状況だろうと速い球を難しいところに打てればいいと考えていました。リスクなんて考えず、とりあえずスマッシュを速く打てばいいと思ってました。難しいショットが決まれば楽しいですからね。でも、それでは強くはなれなくて、リスクコントロールをしていくと安定して勝てるようになるということを大学のサークルで強い先輩方に教わって学んでいきました。こういったリスクコントロール術はあまり教えてくれるところはないし、教本などにもなかなか書いていないことだと思います。さらに、人に言葉で伝えようにもなかなか伝えるのが難しい。そんなわけで数理モデルに落とし込んで普遍的に伝えられるようになればいいのになと思っていて、半年前くらいになんとなく作ってみたというのが、モデル作成の経緯です。

ゲーム理論を学んだことのある方は、「しっぺ返し戦略」や「トリガー戦略」というネーミングに聞き覚えがあったかもしれません。これらは、囚人のジレンマの繰り返しゲームで出てくる戦略です。今回のモデルもゲーム理論ライクなモデルでした。実際のゲーム理論での分析との違いは、ゲーム理論ならば戦略とそれに対応する利得が与えられて、それを最大化するようにプレイヤーが戦略を選択するとして、どのような戦略選択が均衡として落ち着くかを分析するのに対して、今回のモデルは戦略を選ぶとどのような利得(勝ちか負けか)を計算して終わりというものでした。ゲーム理論風に、分析を続けるならば、互いに天邪鬼戦略を選ぶのが均衡になるものだと思われます。

 

LTの時間が5~10分だったのは、発表者にとってはやりやすかったです。超メジャースポーツじゃない場合、スポーツの紹介が必要になってくるので、そこに時間がかかります。自己紹介、スポーツの紹介、分析の紹介とやると、5分じゃ厳しいように思います。研究室や職場、一つのテーマに絞った勉強会のLTであれば、互いに何を知っているかがある程度共有できているので、5分でもいいかもしれませんが、spoanaのように、いろんなスポーツの方がいる場では、5分より10分の方がやりやすいように思えました。

 

 

LT担当があったので、緊張していて、発表全てを集中して聞くことはできていないのですが、以下に、発表中に自分の考えていたこととか絡めながらの感想を書いていきます。

 

まず、全体として、スポーツの分析について、見えるデータの分析と見えない選手のメンタルや思考を炙り出したり、それらに影響を与えたりするのが面白いなって思いました。次に、各々について、振り返りながら記していきます。

spoana.connpass.com

 

エンデュランス競技の話
「自信をもって休むことができる」ってのが良かった。これ大事だなって。こういう非自明なことを明らかにしてプレーを改善していくことが分析の重要な役割の1つなはず。

 

ブラインドサッカーの話
すごいスポーツだと思った。選手らは何が見えているんだろうか。コートの外からの指示者や見えているキーパーはどのように指示を送っているのかが気になった。日本の戦術がカウンター中心から、相手陣地で展開やキーパー軸でのパス回しなどは理に適っているなと思った。ルーズボールの割合のデータを示していたところも良かったです。まだまだ戦略の幅はありそうで、これから盛り上がりそうな競技だと思いました。

 

・相撲の巴戦の話

巴戦の試合回数の理論値とのズレに着目したものでした。
相撲の話といえば、ヤバい経済学でも挙げられていますね。こちらは最終戦の勝ち越し負け越しがかかった試合での勝率のズレを見たものでした。以下がリンクです。

http://pricetheory.uchicago.edu/levitt/Papers/DugganLevitt2002.pdf

この論文は、構成、とくにイントロがとても素晴らしいって話を経済学者の川口先生がツイートしてました。
巴戦の各プレイヤーの勝率については、数学の問題として定式化していて、最初の試合を待つ人が少し不利になるなんて話があって面白いです。

mathtrain.jp

解法1みたいな再帰的な構造がパッと見えるといいなと。解釈としては、一回目待機マンは初戦から絶対に負けられない分だけ損してるんだなと思います。

この確率構造のシミュレーションはspoana運営者の一人のu++さんがやってますね。

upura.hatenablog.com


セーリングのデータ計測の話
自分でデータをスマホで記録するというアイデアが良かったです。たしかにGPSデータや角度のデータなど取得しやすいですし。自分でも動画を撮ったりするけど、自分が扱えるデータに落とし込むのが難しいなと感じていたのでなるほどと。分かりやすいデータにすると、日々の自分の成長が分かって、継続するモチベーションにもなっていいですね。

 

・卓球のTリーグの話
サービス側での得点率とレシーブ側の得点率で、サービスで安定している人の方が強い傾向が見えました。卓球はサービス2回で交代でしたっけ。当然サービスの有利性はスポーツによって違いますが、サービスの重要性は自分がバドミントンをやっていて実感しています。サービスというプレーは、相手の状態によらずルールの範囲内でなら自分の決めた体勢で打てる唯一のショットなんですよ。他のショットは相手の打ったものに依存してしまう。ここを安定させない手はないと思っています。

 

富士登山競争の話
出場資格獲得から難関とあって、ハンター試験か!って心の中でツッコみましたね。つらみ指標の作成は面白い試みだと思いました。高尾山比で、あやしい指標かもしれないということでしたが、高度による空気の補正(METS?)の働きが強すぎたように思えます。そのへんの補正値を調整すると良いかも知れませんね。あと、つらみは瞬発系と持久系で異なるとことかが大事な気もします。個人的には持久系がとてもつらいので。あと、いつか富士山登ってみたいなってなりました、ゆっくりで。


・サポーターの熱量の話
SNSテキストマイニングの感情分析の話。リアルタイムで分析していけるってすごいですね。個人的には会社紹介のところの、知的財産 の裁判における判決の予測が刺さりました。最近、中国の裁判がAIで、なんてニュースを見かけたような気がしたので。


・サッカーのボール保持力の話
スポーツの論文とかあまり知る機会がなかったので、MITの論文とか紹介を受けて、こういうのがあるのかと知れてよかったです。サッカーのようなメジャースポーツはデータや手法が揃っていていいなぁと少し嫉妬したりもしました。

 

・サッカーのポジショナルプレーの話
ボロノイで守備範囲を見たりしていました。トポロジカルデータ解析の話は一度別のところで発表(以下のスライド)を聞いたくらいですが、面白そうな手法だなと気になってはいます。

https://www.slideshare.net/YutakaKuroki/lytics-112543204

ボロノイは商圏の分析とかに使えるかなと思ったりもしました。あと、ボールを持っていないときの動き、オフザボールの動きとか好きなので、そういう分析が出てくるといいな(探せばあるとは思いますが)と。バドミントンをやっていても、ダブルスで球を打っていないペアの動きがコースを制限したりして、実は重要な役割を果たしていると感じています。
少し関係ないですが、ランディ・パウシュという方の最後の授業というタイトルの講演を昔見たときに、こういう考え方を得たのを覚えています。彼が子供時代に、アメフトのコーチにまず教わったのは、アメフトは11人vs11人でやるが、ボールを持っているのは1人だけで、ボールを持っていない21人の動きが実は重要で、その練習をやるんだってことが語られていました。

【全文】ランディ・パウシュ「最後の授業」--余命半年の男が生徒たちに遺した「夢の叶え方」 - ログミーBiz

 

リオ五輪のメダル予想の話
アカデミアの人がこういう場にやってくるのすごくいいですね。ロジスティック回帰でも、汎用的にこれだけの精度が出るんだっていうのがいい。得点比で強さを測るのが肝だったように記憶してます。6月に論文として出版される予定らしいです。

https://search.ieice.org/bin/summary_advpub.php?id=2018EDP7315&category=D&lang=E&abst=

この話とは少し脱線しますが、今回、自分の分析ではやらなかったことの一つとして、点数差に応じて戦略変えるということがありました。というのも、点数差が開いて勝っていれば全力を出さなくてもいいかというものとか、負けている側はもう勝てないからいいかと諦めたりとか、人間の行動が点数差に影響されるだろうなという予測があります。サッカーだと2-0が危ないとかよく言いますよね。個人的な解釈ですが、勝てるもんだと慢心して集中力が欠如しているところで、相手の奮起で得点を入れられると、不安が連鎖したり、各々が負けないよう頑張ろうと頑張るが、統率がズレて穴が開いてさらにやられて、悪い状況が増幅するんだろうなと思っています。バドミントンでも、個人的には、敵のペアの各々が自分より強くても、ペア同士のコミュニケーションがとれていない相手は、サッカーの2-0からの逆転のような状況を作り出すと、大金星がとれたりする気がします。


・部活の話
練習時間とかレギュラーとか人間関係とか難しそうですよね。実は、話を聞きながら北宇治高校吹奏楽部を思い浮かべていました。個人的にも、中高は部活、大学はサークルと所属していたのですが、部活は強くなることという一つの方向をメンバーが向いているがメリットでありデメリットである一方、サークルは各々がいろんな方向を向いているというのがこれまたメリットでありデメリットであるなぁと思いながら過ごしてきました。

 

登壇する機会を与えてくださったSports Analyst Meetup運営の皆さま、参加者や関係者の皆さま、ありがとうございました。今回が初参加でしたが、スポーツの分析の面白さに惹きつけられたので、また自分で分析してみたり、他の方の分析を聴いたりしたいので、また参加していきたいと思います。今後ともよろしくお願いいたします。

 

 

Tokyo.R #76の後日譚々

Tokyo.R #76 @DeNA にてLTを行ってきました。そのスライドがこちらになります。

 

speakerdeck.com

 

発表での伝わりやすさを考慮して、厳密さに関して少し欠落しているところやLDAなのにディリクレ分布の説明が出てこないところ等、至らぬところは目を瞑って、参考文献を読んでいただけると幸いです。

 

コードはこちらの記事にまとめました。

moratoriamuo.hatenablog.com

 

当初のタイトルとしては、『トピックモデル ~江戸時代のレシピを添えて~』というタイトルでした。
人文学オープンデータ共同利用センター

データセット | 人文学オープンデータ共同利用センター
江戸料理レシピデータセットがあり、
これを現代風のレシピにしたものがクックパッドにあるのを発見しました。

cookpad.com


これにトピックモデルを適用してみたが、あまり卵料理が多かったのと、単純に量が少なかった(70件ほど)せいか、あまり面白い結果が出なかったので、急遽路線を現代のレシピに変更して、今回の発表になりました。


それ以前のLT案として、

・『占いは統計学?占い自動生成人工知能 uranAIを作ろう』
データが集められなかったため断念しました。あと、作ってしまうと何かしらまずいことが起きるかもしれないという恐れもありました。
ただ、やってみたい気持ちもあるので、もしも占いのテキストデータを持っている人がいれば、お話ししたいところです。占いではなく、おみくじでも可能です。

・『お酒の失敗談をテキストマイニングする』
少しいじっていたのですが、発表の何かしらの規約に引っかかる可能性を危惧して断念しました。

 

というものがありましたが、これらはボツ案としてお蔵入りです。

 

スライド作成については、Rmarkdownでスライド作ってみたかったのですが、時間の関係上、習得が間に合わず断念し、パワポで作りました。スライドのデザインとしては、クッキング番組風に、「このように作ったものがこちらになります」 みたいなイメージで作りました。色合いもそんな感じで、オレンジ主体になってます。

 

登壇する機会を与えてくださったTokyo.R運営の皆さま、参加者や関係者の皆さま、ありがとうございました。またネタを考えて登壇したいと思っているので、引き続きよろしくお願いいたします。

 

 

トピックモデルで1週間の献立をレコメンドする衝動遊戯:第1ゲーム

 

【ABSTRACT】

料理において、最も面倒くさいことの1つが献立決めである。献立決めのレコメンドエンジンを作れば、この料理への心理的障壁を乗り越えることができる。料理の献立決めのレコメンドエンジンの要件として、同じようなものが続けてレコメンドされないことが求められる。なぜならば、同じものばかりでは飽きてしまうし、栄養の偏りも出てしまうためだ。したがって、この記事では、クックパッドのレシピに、トピックモデルを適用して、レシピをトピックごとにクラスタリングし、クラスタ別で1週間分のレシピを提案するレコメンドエンジンを作成する。

【OBJECTIVE】

料理において、最も面倒くさいことの1つが献立決めである。献立を考える手間を省くため、レコメンドエンジンを作成したい。最も簡単なものとして、レシピ集をランダムに開いて決定する方法が思い浮かぶが、完全にランダムに決定してしまうと、同じものばかり出てしまう可能性が生じる。レシピ集の料理の種類の割合によって出てくるものに偏りが出るのだ。同じものばかりでは飽きてしまうし、栄養のバランスも崩れてしまう。これを避けるために、レシピをクラスタリングして、1週間分の献立を、クラスタが被らないようにしてランダムにレコメンドする方法が考えられる。以降、このようなレコメンドエンジンを作成する。

【METHODS】

f:id:moratoriamuo271:20190305004151p:plain

このような手順でRで作成する。

【RESULTS】

材料確保(スクレイピングする)

使用するライブラリと乱数の種。

library(tidyverse)
library(stringr)
library(rvest)
library(RMeCab)
set.seed(20190302)

クックパッドのカテゴリのリストが
https://cookpad.com/category/list にある。「今日のご飯・おかず」、「お菓子」、「パン」、「離乳食」、「その他」のカテゴリから、「今日のご飯・おかず」を選択する。13214ページあって、レシピ141355品が掲載されている(この件数は日々変動しているので注意)。

14万件全てを取ってくるのは大変なので、50ページ分だけランダムに取ってくることにする。

url <- "https://cookpad.com/category/177?page="
allrepnum <- 1:13214
repset <- sample(allrepnum, 50, replace = FALSE)

ページからレシピの番号(https://cookpad.com/recipe/の後に続く数字)を入手する関数を作成

fget_recipenumlist <- function(url,repset){
  rec <- vector("list", length(repset))
  for(i in 1:length(repset)){
    urlfull <- str_c(url, repset[i] %>% as.character())
    rec[[i]] <- read_html(urlfull) %>% html_nodes(css = "a") %>% 
      html_attr("href") %>% str_subset("/recipe/(\\d\\d\\d\\d)") 
  }
  return(rec)
}

関数を適用する。重複して取られてしまうケースがあるので重複削除。

reclist <- fget_recipenumlist(url, repset)
reclist <- reclist %>% unlist() %>% unique()    #重複を消す

入手したレシピの番号からレシピの中身を入手する関数を作成。ごり押し気味。

fget_cookpads <- function(reclist){
  rec <- vector("list", length(reclist))
  for(i in 1:length(reclist)){
    urlfull <- str_c("https://cookpad.com",reclist[i])
    temp_text  <- read_html(urlfull) %>% html_nodes(css = "p") %>%html_text()
    loc_start <- temp_text %>% str_which("\nメンバー検索") + 1
    loc_end <- temp_text %>% str_which("(\\d\\d)/(\\d\\d)/(\\d\\d)") %>% as.list()
    if(length(loc_end)==0){
      loc_end <- temp_text %>% str_which("お困りの方はこちら") %>% as.list()
    }
    loc_end <- loc_end[[1]] - 1
    rec[[i]] <- temp_text[loc_start:loc_end]
    if(i%%100==0)print(i)    #経過状況を見るため100個おきにプリント
    Sys.sleep(1)
  }
  return(rec)
}

関数を適用する。改行を消すのと、同じレシピは一つの文書にまとめる。

cookpads <- fget_cookpads(reclist)
#改行の\nを消して各レシピを一つの文書にまとめる
cookpads <- purrr::map(cookpads, str_remove_all,"\n") %>% purrr::map(str_c, collapse="")

入手したレシピを各レシピごとで保存する。加えて、レシピを一つの文書にまとめて保存する。

#レシピの番号
reclistnum <- str_extract(reclist, "(\\d){4,}")

#レシピの保存(各レシピ)
for(i in 1:length(cookpads)){
  filename <- str_c("docs_nowrecipe/",reclistnum[i] %>% as.character()) %>% str_c(".txt")
  write(cookpads[[i]],filename)
}

#全レシピ連結で保存
allcookpads <- cookpads %>% unlist() %>% str_c(collapse = "")
write(allcookpads, "allnowdocs.txt")

これにより、500件のレシピが保存できた。

材料確認(ワードクラウドで遊んでみる)

library(tidyverse)
library(stringr)
library(RMeCab)
library(tm)
library(topicmodels)

どんな単語があるかを確認。自立語(動詞や形容詞)で見るケースと名詞で見るケース。

all_freq <- RMeCabFreq("allnowdocs.txt")
## file = allnowdocs.txt 
## length = 4110
all_freq %>% filter(Info2=="自立") %>% arrange(-Freq) %>% head(20)
##          Term  Info1 Info2 Freq
## 1        する   動詞  自立 1366
## 2      入れる   動詞  自立  897
## 3        切る   動詞  自立  412
## 4      加える   動詞  自立  332
## 5      炒める   動詞  自立  325
## 6      混ぜる   動詞  自立  285
## 7        なる   動詞  自立  193
## 8      かける   動詞  自立  174
## 9  出来上がる   動詞  自立  171
## 10       焼く   動詞  自立  167
## 11     茹でる   動詞  自立  153
## 12     食べる   動詞  自立  138
## 13       熱す   動詞  自立  117
## 14       煮る   動詞  自立  116
## 15       軽い 形容詞  自立  113
## 16     大きい 形容詞  自立  107
## 17       作る   動詞  自立   98
## 18       取る   動詞  自立   87
## 19       盛る   動詞  自立   81
## 20     つける   動詞  自立   80
all_freq %>% filter(Info1=="名詞"&Info2!="数"&Info2!="接尾") %>% arrange(-Freq) %>% head(20)
##          Term Info1    Info2 Freq
## 1  フライパン  名詞     一般  307
## 2          塩  名詞     一般  272
## 3          火  名詞     一般  254
## 4          水  名詞     一般  241
## 5          油  名詞     一般  182
## 6          ♪  名詞 サ変接続  178
## 7          鍋  名詞     一般  163
## 8          肉  名詞     一般  160
## 9          卵  名詞     一般  136
## 10          +  名詞 サ変接続  131
## 11          U  名詞     一般  130
## 12     玉ねぎ  名詞     一般  129
## 13       よう  名詞   非自立  129
## 14         皮  名詞     一般  121
## 15     レシピ  名詞     一般  117
## 16         味  名詞     一般  117
## 17       調味  名詞 サ変接続  114
## 18       材料  名詞     一般  113
## 19         皿  名詞     一般  112
## 20       好み  名詞     一般  106

名詞と自立語でワードクラウドを作ってみる。(警告が現れるが、楽しく可視化したいだけなので無視している。)

set.seed(271)
corpus <- docDF("alldocs.txt", type = 1)
## file_name =  ./alldocs.txt opened
## number of extracted terms = 1166
## now making a data frame. wait a while!
corpus <- corpus[81:nrow(corpus),]    #記号とか余計なものを消した
options(warn=-1)    #ワードクラウドの警告が邪魔なので消した
corpus1 <- corpus %>% filter(POS1 == "名詞"&POS2=="一般"| POS2 == "自立") 
wordcloud::wordcloud(corpus1$TERM, corpus1$alldocs.txt, min.freq= 5, scale=c(15,0.5),
                     family ="JP1", random.color=TRUE, color = rainbow(5))

options(warn=0)

名詞でワードクラウドを作ってみる。

options(warn=-1)
corpus2 <- corpus %>% filter(POS1 == "名詞"&POS2=="一般")
wordcloud::wordcloud(corpus2$TERM, corpus2$alldocs.txt, min.freq= 5, scale=c(15,0.5),
                     family ="JP1", random.color=TRUE, color = rainbow(5))

options(warn=0)

調理(トピックモデル)

トピックモデルの1つのLDA(Latent Dirichlet Allocation)を用いてクラスタリングを行う。

名詞と動詞に絞って、文書ターム行列を作成。記号等も削除している。削除の仕方はごり押し気味。

dir_recipe <- "docs_nowrecipe"
DTM <- docDF(dir_recipe, type = 1)
DTM <- DTM[512:4106,]

LDAを適用するために文書ターム行列の形に整形する。

DTM_nouns_verbs <- DTM %>% filter( (POS1 == "名詞"&POS2 == "一般") | POS2 == "動詞" )%>% 
  select(-c(POS1,POS2)) %>% 
  data.frame(row.names = 1) %>% t() %>% 
  as.DocumentTermMatrix(weighting = weightTf)

後でperplexity計算するため8:2でtrainとtestに分割する。

#500*1720の行列
DTM_nouns_verbs_train <- DTM_nouns_verbs[1:400,]
DTM_nouns_verbs_test <- DTM_nouns_verbs[401:500,]

モデル選択(トピック数の決定)の方法として

  • パープレキシティ perplexityによる評価 負の対数尤度から計算される値。testデータを使用。低い方が良いモデル。

  • {ldatuning}パッケージの使用 4つの論文で提案された指標で評価する。

  • 変分下限でモデル評価 変分ベイズを使っているならば評価に使える。

  • ディリクレ過程でトピック数もモデルに組み込む 階層ディリクレ過程を用いるとトピック数の推定が可能。勉強不足のため、これについてはまだよく分かっていない。pythonのgenismというものを使うとできるかもしれない。以下のブログが参考になるかもしれない。 https://qiita.com/yamano357/items/90863cf15326ebd19231

以下では、{ldatuning}パッケージで、だいたいのあたりをつけて、その中からperplexityで評価して、トピック数の選択を行う。
(トピック数のベストな設定方法は結局のところ分かっていないので、もう少し調べていきたいところ。混合ユニグラムモデルまでは過去にシミュレーションデータが作れているので、LDAのシミュレーションデータを作って各指標をチェックする記事をそのうち書きたいと思っている。)

ここでは、{ldatuning}のパッケージでトピック数を2から97のトピック数を5個刻みにして評価する。

result <- ldatuning::FindTopicsNumber(
  DTM_nouns_verbs_train,
  topics = seq(from = 2, to = 97, by = 5),
  metrics = c("Griffiths2004", "CaoJuan2009", "Arun2010", "Deveaud2014"),
  method = "Gibbs",
  control = list(seed = 777),
  mc.cores = 4L,
  verbose = TRUE
)
## fit models... done.
## calculate metrics:
##   Griffiths2004... done.
##   CaoJuan2009... done.
##   Arun2010... done.
##   Deveaud2014... done.
ldatuning::FindTopicsNumber_plot(values = result)