【R/Stan】潜在変数モデリングで潜在的に評価の高い南アジア料理の店を探す
突然ですが、最近のマイブームは都内の美味しい南アジア料理屋を探し回ることです。とはいっても、かなりの店舗さんが存在しますのでスパッと「ここいいかも～」みたいなお店を探し当てるのは中々苦労します。その中でも役に立つのが皆さんご贔屓の「食べログ」ですが、その評価スコアについて「点数の基準がようわからん」とか「何か点数にカラクリがあるんじゃね」のような疑念が絶えず向けられており、私自身も常に疑いながらdigってます。
そこで、いわゆる食べログスコアだけを参考にするだけではなくて「口コミ数」「ブックマーク数」「食事代」「立地」などの情報を利用することで潜在的な評価点を推定すれば、良いお店を探し当てることができるのでは？！とふと思い立ちました。ここで皆さんのお待ちかね、潜在変数モデルの出番ですよね。というわけで、今回もRとStanをこねくり回して潜在的に評価の高い南アジア料理屋さんはどこなのかを掘り当てていきます。
ちなみに、食べログのデータはサイトからスクレイプして持ってきたのですが、諸々説明してると長くなるので後程記事にしようと思います。気になる方はそちらをお読みください。今回はその辺端折ります、スミマセン。結果だけ知りたい人は下記の目次から飛んでください。
準備
まずは必要なパッケージとデータを読み込みます。スクレイプで取得したデータはsouthasia_data.RDataとして保存しました。
library(tidyverse) # データ加工 library(cmdstanr) # Stanを動かす library(posterior) # cmdstanの推定結果をハンドリングするやつ load("southasia_data.RData") # 南アジア料理屋の食べログ情報が入ったデータ
最近になってrstanからcmdstanrに乗り換えました。というのも、再コンパイルの必要がなかったりセッションがクラッシュすることが少なくなったりと色々なメリットがあったのでそうしました。また、WSL2でRstudio serverを動かすようになったので最近かなり快適なStanライフが送れてる気がします。
データの確認をしてみます。
dim(southasia_data)
[1] 1998 8
店舗数は全部で1998店舗といったところでしょうか。
head(southasia_data)
# A tibble: 6 × 8 name city area score nreview bookmark lunch_bud dinner_bud <chr> <chr> <chr> <dbl> <dbl> <dbl> <chr> <chr> 1 Kalpasi 世田谷区 千歳船橋 3.9 165 30774 ￥2,000～￥2,999 ￥3,000～￥3,999 2 シバカリーワラ 世田谷区 三軒茶屋 3.87 475 52270 ￥1,000～￥1,999 ￥1,000～￥1,999 3 カッチャル バッチャル 豊島区 新大塚 3.85 641 67299 ￥1,000～￥1,999 ￥3,000～￥3,999 4 砂の岬 世田谷区 桜新町 3.84 239 30972 ￥1,000～￥1,999 NA 5 やっぱりインディア 豊島区 大塚 3.83 415 24542 ￥1,000～￥1,999 ￥2,000～￥2,999 6 スパイスバル コザブロ 文京区 本駒込 3.83 159 15341 ￥1,000～￥1,999 ￥4,000～￥4,999
変数は以下の通りです。
name- 店舗名city- 市区町村area- 最寄り駅score- 食べログスコアnreview- 口コミ数bookmark- ブックマーク数lunch_bud,dinner_bud- 食事代 (ランチ、ディナー)
とりあえず、市区町村ごとにどれ程の店舗があるのかを確認してみます。せっかくなので東京都の地図上にプロットしてみましょう。{jpndistrict}パッケージ を利用すれば東京都のsfデータを取得することができます。
# 東京都の地図データ tokyo <- jpndistrict::jpn_pref(13) # 島嶼部を除外 tokyo <- tokyo[-(54:nrow(tokyo)), ] # 市区町村ごとに店舗数をカウント count_shop <- southasia_data %>% group_by(city) %>% summarise(n = n()) # プロット tokyo %>% left_join(count_shop, by = "city") %>% mutate(n = if_else(is.na(n), as.integer(0), n)) %>% ggplot() + geom_sf(aes(fill = n)) + theme_bw() + scale_fill_gradient(low = "white", high = "springgreen4", name = "Count") + ggtitle("Number of restaurants per municipality")
都心に相当集中してますが、意外と八王子にも郊外にしては店舗数がありますね。
生成モデル
潜在的評価をいかにして推定するかについて、以下の生成モデルを考えます。
\begin{align}
Y_{ij} \sim \mathcal{N}(\alpha_j + \beta_j \theta_i, \kappa_\mu).
\end{align}
$i = 1, \dots, I$とし、 料理店を表すインデックスとします。また、$j = 1, \dots ,3$は三つの評価項目を表し、「食べログスコア」「口コミ数」「ブックマーク数」です。$Y_{ij}$ は料理店$i$の評価項目$j$での獲得ポイントとします。
上式は、ほとんど項目反応理論 (IRT) の式形と同じです。$\alpha_j$は切片項で$\beta_j$は因子負荷量だとします。また、$\theta_i$がお店の潜在的評価です。IRT風に言えば、それぞれ「困難度」「識別度」「能力」といったところです。$\alpha_j$は項目のポイントを取りやすさを表すパラメータであり、$\beta_j$は潜在的評価$\theta_i$が高い店と低い店を切り分けるような項目のパラメータです。
ここで、項目のパラメータである$\alpha_j$, $\beta_j$は階層構造があることを仮定します。$c$を市区町村、$l$をランチの食事代カテゴリー、$d$をディナーの食事代カテゴリーとします。これらを導入することによって、評価項目のパラメータが市区町村やランチ・ディナー食事代カテゴリーごとに別々に生成されるようにします。さらに、市区町村ごとや食事代ごとによる項目評価の重みの違いなどを考慮したモデルになってます。具体的に書くと、
\begin{align}
\alpha_j &= \delta_{j} + \alpha_{cj} + \alpha_{lj} + \alpha_{dj},\\
\beta_j &= \gamma_{j} + \beta_{cj} + \beta_{lj} + \beta_{dj}.
\end{align}
となります。$\delta_{j}$ は評価項目 $j$ の全体平均項目難易度、$\alpha_{cj}$ は市区町村 $c$ における項目難易度、$\alpha_{lj}$ はランチ代カテゴリー$l$における項目難易度、$\alpha_{dj}$はディナー代カテゴリー$d$における項目難易度です。同様にして、$\beta$に関してもそれぞれの層ごとの識別度になります。
また、階層的な事前分布を以下で与えます。
\begin{align*}
\delta_j &\sim \mathcal{N}(0, 10), \ \ \ \gamma_j \sim \mathcal{N}(0, 10), \\
\alpha_{cj} &\sim \mathcal{N} (0, \kappa_{\alpha_C}), \ \ \ \beta_{cj} \sim \mathcal{N}(0, \kappa_{\beta_C}),\\
\alpha_{lj} &\sim \mathcal{N}(0, \kappa_{\alpha_L}), \ \ \ \beta_{lj} \sim \mathcal{N}(0, \kappa_{\beta_L}),\\
\alpha_{dj} &\sim \mathcal{N}(0, \kappa_{\alpha_D}),\ \ \ \beta_{dj} \sim \mathcal{N}(0, \kappa_{\beta_D}).
\end{align*}
$\kappa_{\alpha}, \kappa_{\beta}$はそれぞれの層ごとの標準偏差です。それぞれの層ごとに評価項目パラメータを生成します。また、潜在的評価$\theta_i$の事前分布を
\begin{equation*}
\theta_i \sim \mathcal{N}(0, 1).
\end{equation*}
とし、$\kappa_\mu$の事前分布を
\begin{equation*}
\kappa_\mu \sim Ga(1, 1).
\end{equation*}
とします。また、食べログスコア $S$と比較しやすくするため、$[min(S), max(S)]$の範囲で$\theta_i$を正規化します。
\begin{equation*}
\theta_i^\star = \frac{(\theta_i - min(\theta)) (max(S) - min(S)) }{max(\theta) - min(\theta)} + min(S).
\end{equation*}
これらをStanファイル上で記述してみます (tabelog.stan)。
data {
int N; //# of observations
int I; //# of shops
int J; //# of evaluation items
int C; //# of cities
int L; //# of lunch budget categories
int D; //# of dinner budget categories
int i[N]; //shop index
int j[N]; //evaluation item index
int c[N]; //city index
int l[N]; //lunch budget category index
int d[N]; //dinner budget category index
vector[N] Y; //evaluations
real max_tabelog; //maximum tabelog score for standardizing
real min_tabelog; //minimum tabelog score for standardizing
}
parameters {
vector[I] theta; //latent value
vector[J] delta; //mean difficulty
vector[J] alpha_C[C]; //city level difficulty
vector[J] alpha_L[L]; //lunch budget category level difficulty
vector[J] alpha_D[D]; //dinner budget category level difficulty
vector[J] gamma; //mean discrimination
vector[J] beta_C[C]; //city level discrimination
vector[J] beta_L[L]; //lunch budget discrimination
vector[J] beta_D[D]; //dinner budget discrimination
real<lower=0> sigma_aC; //sd for city level alpha
real<lower=0> sigma_aL; //sd for lunch budget level alpha
real<lower=0> sigma_aD; //sd for dinner budget level alpha
real<lower=0> sigma_bC; //sd for city level beta
real<lower=0> sigma_bL; //sd for lunch budget level beta
real<lower=0> sigma_bD; //sd for dinner budget level beta
real<lower=0> sigma_mu; //sd for evaluation score
}
transformed parameters {
vector[J] alpha; //difficulty
vector[J] beta; //discrinimation
vector[N] mu; //linear predictor
for (n in 1:N) {
alpha[j[n]] = delta[j[n]] + alpha_C[c[n], j[n]] + alpha_L[l[n], j[n]] + alpha_D[d[n], j[n]];
beta[j[n]] = gamma[j[n]] + beta_C[c[n], j[n]] + beta_L[l[n], j[n]] + beta_D[d[n], j[n]];
if (beta[j[n]] < 0) {
beta[j[n]] = beta[j[n]] * -1;
}
mu[n] = alpha[j[n]] + beta[j[n]] * theta[i[n]];
}
}
model {
//multilevel priors
for (cc in 1:C) {
alpha_C[cc] ~ normal(0, sigma_aC);
beta_C[cc] ~ normal(0, sigma_bC);
}
for (ll in 1:L) {
alpha_L[ll] ~ normal(0, sigma_aL);
beta_L[ll] ~ normal(0, sigma_bL);
}
for (dd in 1:D) {
alpha_D[dd] ~ normal(0, sigma_aD);
beta_D[dd] ~ normal(0, sigma_bD);
}
//priors
theta ~ std_normal();
delta ~ normal(0, 10);
gamma ~ normal(0, 10);
sigma_mu ~ gamma(1, 1);
//likelihood
Y ~ normal(mu, sigma_mu);
}
generated quantities {
//standardization
vector[I] theta_std;
theta_std = (theta - min(theta)) / (max(theta) - min(theta)) * (max_tabelog - min_tabelog) + min_tabelog;
}
Stanを動かす
それでは早速Stanによる推定に移っていこうと思います。ここで、Stanに送るためのデータが必要になってくるのでスパッと作っていきましょう。
正直必要のない作業だとは思いますが行列化した方が色々やりやすいので、とりあえずは先ほどのデータフレームを行列に変換します。また、口コミ数とブックマーク数はバラツキが非常に大きいので自然対数変換します。
mat <- southasia_data %>% mutate(ln_nreview = log(nreview), ln_bookmark = log(bookmark)) %>% select(name, score, ln_nreview, ln_bookmark) %>% as.matrix() sname <- mat[, 1] # 店舗名を抜き出す mat <- mat[, -1] # 一列目が店舗名列なので削除 mat <- apply(mat, 2, as.numeric) # なぜかcharacter matrix扱いになっているので全部数値に変換します rownames(mat) <- sname # 店舗名で行をラベリング
必要な変数を作り、リストにまとめます。
Y <- as.numeric(mat) # スコア、口コミ数、ブックマーク数をひとまとめ N <- length(Y) # 観測数 i <- rep(1:nrow(mat), times = ncol(mat)) # 店舗インデックス I <- max(i) # 店舗数 j <- rep(1:ncol(mat), each = nrow(mat)) # 項目インデックス J <- max(j) # 項目数 c <- rep(southasia_data$city_id, times = J) # 市区町村インデックス C <- max(c) # 市区町村数 l <- rep(southasia_data$lunch_bud_id, times = J) # ランチ代カテゴリーインデックス L <- max(l) #ランチ代カテゴリー数 d <- rep(southasia_data$lunch_bud_id, times = J) # ディナー代カテゴリーインデックス D <- max(d) # ディナー代カテゴリー数 max_tabelog <- max(southasia_data$score) # 食べログ最高スコア（標準化のため） min_tabelog <- min(southasia_data$score) # 食べログ最低スコア（同上） data <- list(N = N, I = I, J = J, C = C, L = L, D = D, i = i, j = j, c = c, l = l, d = d, Y = Y, max_tabelog = max_tabelog, min_tabelog = min_tabelog)
次に、上記で作成したtabelog.stanをコンパイルします。
model <- cmdstan_model("tablelog.stan", cpp_options = list( stan_threads = TRUE # スレッド数を指定 ))
読み込めたので分析を実行していこうと思います。普通のサンプリングだと計算がかなり遅いので、今回は自動変分ベイズ (ADVI) で計算スピードを爆上げします。モデル名$variationalで{cmdstanr}でADVIを実装できます。
fit <- model$variational(data, seed = 1, threads = parallel::detectCores())
推定には10秒もかからなかったと思います。結果を抽出していきます。{cmdstanr}ではfit$summary("パラメータ")で任意の推定値を抜き出すことができますが、デフォルトで90%区間しか取ってきてくれないので95%が欲しいときは自力で計算する必要があります。fit$draws("パラメータ")で近似分布からのパラメータ推定値 1000個 (デフォルト) を抜き出してくれます。その推定値を基にして、95%信用区間を計算していきます。まず、下側・上側を計算する関数を定義します。
lower <- function(x) quantile(x, probs = .025) upper <- function(x) quantile(x, probs = .975)
潜在的な評価である$\theta_i$を$[min(S), max(S)]$範囲で正規化したパラメータである$\theta_i^\star$の推定値を抜き出してみましょう。
theta <- fit$draws("theta_std") theta_est <- tibble(name = sname, station = southasia_data$area, city = southasia_data$city, lunch_bud = southasia_data$lunch_bud, dinner_bud = southasia_data$dinner_bud, lwr = apply(theta, 2, lower), med = apply(theta, 2, median), upr = apply(theta, 2, upper), tabelog_score = southasia_data$score)
結果
とりあえず推定指標と食べログスコアの相関をプロットしてみます。
cor <- cor.test(theta_est$med, theta_est$tabelog_score) theta_est %>% ggplot(aes(x = med, y = tabelog_score)) + geom_point() + theme_bw() + xlab("Estimated latent score") + ylab("Tabelog score") + annotate(geom = "text", label = paste("italic(r) == ", round(cor$estimate, 3)), x = 3.75, y = 3.2, parse = TRUE)
まあ、特に何か得られるわけではないですがある程度の相関はありそうです。それでは推定から得られた潜在的に評価が高い南アジア料理店トップ20を以下に発表していこうと思います。
| 順位 | 店舗名 | 最寄り駅 | 潜在的評価 | ランチ代 | ディナー代 |
|---|---|---|---|---|---|
| 1 | コチン ニヴァース | 西新宿五丁目 | 3.900 | NA | ￥1,000～￥1,999 |
| 2 | ケララバワン | 練馬 | 3.869 | ￥1,000～￥1,999 | ￥2,000～￥2,999 |
| 3 | シャングリーラ 蒲田店 | 蒲田 | 3.851 | ～￥999 | ￥1,000～￥1,999 |
| 4 | カッチャル バッチャル | 新大塚 | 3.849 | ￥1,000～￥1,999 | ￥3,000～￥3,999 |
| 5 | タンドールバル カマルプール 木場店 | 木場 | 3.847 | ￥1,000～￥1,999 | NA |
| 6 | 南印度ダイニング ポンディバワン 武蔵新田 | 武蔵新田 | 3.846 | ～￥999 | ￥2,000～￥2,999 |
| 7 | ハリマ・ケバブ・ビリヤニ | 稲荷町 | 3.830 | ～￥999 | ￥2,000～￥2,999 |
| 8 | スパイスカフェ | 押上 | 3.817 | ￥1,000～￥1,999 | ￥6,000～￥7,999 |
| 9 | アーンドラ・キッチン | 御徒町 | 3.815 | ～￥999 | ￥2,000～￥2,999 |
| 10 | クンビラ | 恵比寿 | 3.808 | ～￥999 | ￥4,000～￥4,999 |
| 11 | シバカリーワラ | 三軒茶屋 | 3.806 | ￥1,000～￥1,999 | ￥1,000～￥1,999 |
| 12 | サルマ ティッカアンドビリヤニ | 品川 | 3.791 | ～￥999 | ￥2,000～￥2,999 |
| 13 | カリーライス専門店エチオピア 本店 | 神保町 | 3.789 | ￥1,000～￥1,999 | ～￥999 |
| 14 | カーン・ケバブ・ビリヤニ | 新橋 | 3.786 | ～￥999 | ￥2,000～￥2,999 |
| 15 | ネパール民族料理 アーガン | 新大久保 | 3.785 | ～￥999 | ￥2,000～￥2,999 |
| 16 | ダバ インディア | 京橋 | 3.780 | ￥1,000～￥1,999 | ￥2,000～￥2,999 |
| 17 | ヴェヌス サウス インディアン ダイニング 錦糸町店 | 錦糸町 | 3.778 | ￥1,000～￥1,999 | ￥1,000～￥1,999 |
| 18 | タンブリン カレー&バー | 北千住 | 3.777 | ￥1,000～￥1,999 | ￥2,000～￥2,999 |
| 19 | Sajilo Cafe | 吉祥寺 | 3.775 | ￥1,000～￥1,999 | ￥1,000～￥1,999 |
| 20 | アーンドラ・ダイニング 銀座 | 銀座一丁目 | 3.773 | ～￥999 | ￥2,000～￥2,999 |
うーん。妥当な順位なのかわからんな～。実際に推定された指標の妥当性を確かめる方法は実際に食べるしかないので、このリストに上がってきたってことで行ってみようっていう気にはなりました。
結論を言うと、何もわかりませんでした！
おわりに
生成モデルを作って都内南アジア料理店の潜在的評価を推定してみましたが、結局数値だけ見ても何もわからないので実際に足を運んで確かめるべしという素晴らしい知見を得ることができました。
ちなみに、私の推しは清澄白河にあるナンディニさんと西早稲田のアプサラさんです。前者は平日のミールスセットが素晴らしくて後者はバナナリーフのスリランカカレーが最高においしいです。最近ミールスに異常にハマっているので良いお店があれば知らせてください。以上。
